1 2 #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 3 #include <petscdmshell.h> 4 #include <petscdmda.h> 5 #include <petscviewer.h> 6 #include <petscdraw.h> 7 8 /* Logging support */ 9 PetscClassId TS_CLASSID, DMTS_CLASSID; 10 PetscLogEvent TS_Step, TS_PseudoComputeTimeStep, TS_FunctionEval, TS_JacobianEval; 11 12 const char *const TSExactFinalTimeOptions[] = {"STEPOVER","INTERPOLATE","MATCHSTEP","TSExactFinalTimeOption","TS_EXACTFINALTIME_",0}; 13 14 struct _n_TSMonitorDrawCtx { 15 PetscViewer viewer; 16 PetscDrawAxis axis; 17 Vec initialsolution; 18 PetscBool showinitial; 19 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */ 20 PetscBool showtimestepandtime; 21 int color; 22 }; 23 24 #undef __FUNCT__ 25 #define __FUNCT__ "TSSetFromOptions" 26 /*@ 27 TSSetFromOptions - Sets various TS parameters from user options. 28 29 Collective on TS 30 31 Input Parameter: 32 . ts - the TS context obtained from TSCreate() 33 34 Options Database Keys: 35 + -ts_type <type> - TSEULER, TSBEULER, TSSUNDIALS, TSPSEUDO, TSCN, TSRK, TSTHETA, TSGL, TSSSP 36 . -ts_save_trajectory - checkpoint the solution at each time-step 37 . -ts_max_steps <maxsteps> - maximum number of time-steps to take 38 . -ts_final_time <time> - maximum time to compute to 39 . -ts_dt <dt> - initial time step 40 . -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 41 . -ts_max_snes_failures <maxfailures> - Maximum number of nonlinear solve failures allowed 42 . -ts_max_reject <maxrejects> - Maximum number of step rejections before step fails 43 . -ts_error_if_step_fails <true,false> - Error if no step succeeds 44 . -ts_rtol <rtol> - relative tolerance for local truncation error 45 . -ts_atol <atol> Absolute tolerance for local truncation error 46 . -ts_adjoint_solve <yes,no> After solving the ODE/DAE solve the adjoint problem (requires -ts_save_trajectory) 47 . -ts_monitor - print information at each timestep 48 . -ts_monitor_lg_timestep - Monitor timestep size graphically 49 . -ts_monitor_lg_solution - Monitor solution graphically 50 . -ts_monitor_lg_error - Monitor error graphically 51 . -ts_monitor_lg_snes_iterations - Monitor number nonlinear iterations for each timestep graphically 52 . -ts_monitor_lg_ksp_iterations - Monitor number nonlinear iterations for each timestep graphically 53 . -ts_monitor_sp_eig - Monitor eigenvalues of linearized operator graphically 54 . -ts_monitor_draw_solution - Monitor solution graphically 55 . -ts_monitor_draw_solution_phase <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom 56 . -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction() 57 . -ts_monitor_solution_binary <filename> - Save each solution to a binary file 58 . -ts_monitor_solution_vtk <filename.vts> - Save each time step to a binary file, use filename-%%03D.vts 59 - -ts_monitor_envelope - determine maximum and minimum value of each component of the solution over the solution time 60 61 Developer Note: We should unify all the -ts_monitor options in the way that -xxx_view has been unified 62 63 Level: beginner 64 65 .keywords: TS, timestep, set, options, database 66 67 .seealso: TSGetType() 68 @*/ 69 PetscErrorCode TSSetFromOptions(TS ts) 70 { 71 PetscBool opt,flg,tflg; 72 PetscErrorCode ierr; 73 PetscViewer monviewer; 74 char monfilename[PETSC_MAX_PATH_LEN]; 75 SNES snes; 76 TSAdapt adapt; 77 PetscReal time_step; 78 TSExactFinalTimeOption eftopt; 79 char dir[16]; 80 const char *defaultType; 81 char typeName[256]; 82 83 PetscFunctionBegin; 84 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 85 ierr = PetscObjectOptionsBegin((PetscObject)ts);CHKERRQ(ierr); 86 if (((PetscObject)ts)->type_name) defaultType = ((PetscObject)ts)->type_name; 87 else defaultType = TSEULER; 88 89 ierr = TSRegisterAll();CHKERRQ(ierr); 90 ierr = PetscOptionsFList("-ts_type", "TS method"," TSSetType", TSList, defaultType, typeName, 256, &opt);CHKERRQ(ierr); 91 if (opt) { 92 ierr = TSSetType(ts, typeName);CHKERRQ(ierr); 93 } else { 94 ierr = TSSetType(ts, defaultType);CHKERRQ(ierr); 95 } 96 97 /* Handle generic TS options */ 98 if (ts->trajectory) tflg = PETSC_TRUE; 99 else tflg = PETSC_FALSE; 100 ierr = PetscOptionsBool("-ts_save_trajectory","Save the solution at each timestep","TSSetSaveTrajectory",tflg,&tflg,NULL);CHKERRQ(ierr); 101 if (tflg) {ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr);} 102 if (ts->adjoint_solve) tflg = PETSC_TRUE; 103 else tflg = PETSC_FALSE; 104 ierr = PetscOptionsBool("-ts_adjoint_solve","Solve the adjoint problem immediately after solving the forward problem","",tflg,&tflg,&flg);CHKERRQ(ierr); 105 if (flg) { 106 if (!ts->trajectory) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONGSTATE,"Must also set -ts_save_trajectory when setting -ts_adjoint_solve"); 107 ts->adjoint_solve = tflg; 108 ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr); 109 } 110 ierr = PetscOptionsInt("-ts_max_steps","Maximum number of time steps","TSSetDuration",ts->max_steps,&ts->max_steps,NULL);CHKERRQ(ierr); 111 ierr = PetscOptionsReal("-ts_final_time","Time to run to","TSSetDuration",ts->max_time,&ts->max_time,NULL);CHKERRQ(ierr); 112 ierr = PetscOptionsReal("-ts_init_time","Initial time","TSSetTime",ts->ptime,&ts->ptime,NULL);CHKERRQ(ierr); 113 ierr = PetscOptionsReal("-ts_dt","Initial time step","TSSetTimeStep",ts->time_step,&time_step,&flg);CHKERRQ(ierr); 114 if (flg) { 115 ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr); 116 } 117 ierr = PetscOptionsEnum("-ts_exact_final_time","Option for handling of final time step","TSSetExactFinalTime",TSExactFinalTimeOptions,(PetscEnum)ts->exact_final_time,(PetscEnum*)&eftopt,&flg);CHKERRQ(ierr); 118 if (flg) {ierr = TSSetExactFinalTime(ts,eftopt);CHKERRQ(ierr);} 119 ierr = PetscOptionsInt("-ts_max_snes_failures","Maximum number of nonlinear solve failures","TSSetMaxSNESFailures",ts->max_snes_failures,&ts->max_snes_failures,NULL);CHKERRQ(ierr); 120 ierr = PetscOptionsInt("-ts_max_reject","Maximum number of step rejections before step fails","TSSetMaxStepRejections",ts->max_reject,&ts->max_reject,NULL);CHKERRQ(ierr); 121 ierr = PetscOptionsBool("-ts_error_if_step_fails","Error if no step succeeds","TSSetErrorIfStepFails",ts->errorifstepfailed,&ts->errorifstepfailed,NULL);CHKERRQ(ierr); 122 ierr = PetscOptionsReal("-ts_rtol","Relative tolerance for local truncation error","TSSetTolerances",ts->rtol,&ts->rtol,NULL);CHKERRQ(ierr); 123 ierr = PetscOptionsReal("-ts_atol","Absolute tolerance for local truncation error","TSSetTolerances",ts->atol,&ts->atol,NULL);CHKERRQ(ierr); 124 125 #if defined(PETSC_HAVE_SAWS) 126 { 127 PetscBool set; 128 flg = PETSC_FALSE; 129 ierr = PetscOptionsBool("-ts_saws_block","Block for SAWs memory snooper at end of TSSolve","PetscObjectSAWsBlock",((PetscObject)ts)->amspublishblock,&flg,&set);CHKERRQ(ierr); 130 if (set) { 131 ierr = PetscObjectSAWsSetBlock((PetscObject)ts,flg);CHKERRQ(ierr); 132 } 133 } 134 #endif 135 136 /* Monitor options */ 137 ierr = PetscOptionsString("-ts_monitor","Monitor timestep size","TSMonitorDefault","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 138 if (flg) { 139 ierr = PetscViewerASCIIOpen(PetscObjectComm((PetscObject)ts),monfilename,&monviewer);CHKERRQ(ierr); 140 ierr = TSMonitorSet(ts,TSMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 141 } 142 ierr = PetscOptionsString("-ts_monitor_python","Use Python function","TSMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 143 if (flg) {ierr = PetscPythonMonitorSet((PetscObject)ts,monfilename);CHKERRQ(ierr);} 144 145 ierr = PetscOptionsName("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",&opt);CHKERRQ(ierr); 146 if (opt) { 147 TSMonitorLGCtx ctx; 148 PetscInt howoften = 1; 149 150 ierr = PetscOptionsInt("-ts_monitor_lg_timestep","Monitor timestep size graphically","TSMonitorLGTimeStep",howoften,&howoften,NULL);CHKERRQ(ierr); 151 ierr = TSMonitorLGCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);CHKERRQ(ierr); 152 ierr = TSMonitorSet(ts,TSMonitorLGTimeStep,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);CHKERRQ(ierr); 153 } 154 ierr = PetscOptionsName("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",&opt);CHKERRQ(ierr); 155 if (opt) { 156 TSMonitorLGCtx ctx; 157 PetscInt howoften = 1; 158 159 ierr = PetscOptionsInt("-ts_monitor_lg_solution","Monitor solution graphically","TSMonitorLGSolution",howoften,&howoften,NULL);CHKERRQ(ierr); 160 ierr = TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&ctx);CHKERRQ(ierr); 161 ierr = TSMonitorSet(ts,TSMonitorLGSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);CHKERRQ(ierr); 162 } 163 ierr = PetscOptionsName("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",&opt);CHKERRQ(ierr); 164 if (opt) { 165 TSMonitorLGCtx ctx; 166 PetscInt howoften = 1; 167 168 ierr = PetscOptionsInt("-ts_monitor_lg_error","Monitor error graphically","TSMonitorLGError",howoften,&howoften,NULL);CHKERRQ(ierr); 169 ierr = TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&ctx);CHKERRQ(ierr); 170 ierr = TSMonitorSet(ts,TSMonitorLGError,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);CHKERRQ(ierr); 171 } 172 ierr = PetscOptionsName("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",&opt);CHKERRQ(ierr); 173 if (opt) { 174 TSMonitorLGCtx ctx; 175 PetscInt howoften = 1; 176 177 ierr = PetscOptionsInt("-ts_monitor_lg_snes_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGSNESIterations",howoften,&howoften,NULL);CHKERRQ(ierr); 178 ierr = TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);CHKERRQ(ierr); 179 ierr = TSMonitorSet(ts,TSMonitorLGSNESIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);CHKERRQ(ierr); 180 } 181 ierr = PetscOptionsName("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",&opt);CHKERRQ(ierr); 182 if (opt) { 183 TSMonitorLGCtx ctx; 184 PetscInt howoften = 1; 185 186 ierr = PetscOptionsInt("-ts_monitor_lg_ksp_iterations","Monitor number nonlinear iterations for each timestep graphically","TSMonitorLGKSPIterations",howoften,&howoften,NULL);CHKERRQ(ierr); 187 ierr = TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,300,300,howoften,&ctx);CHKERRQ(ierr); 188 ierr = TSMonitorSet(ts,TSMonitorLGKSPIterations,ctx,(PetscErrorCode (*)(void**))TSMonitorLGCtxDestroy);CHKERRQ(ierr); 189 } 190 ierr = PetscOptionsName("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",&opt);CHKERRQ(ierr); 191 if (opt) { 192 TSMonitorSPEigCtx ctx; 193 PetscInt howoften = 1; 194 195 ierr = PetscOptionsInt("-ts_monitor_sp_eig","Monitor eigenvalues of linearized operator graphically","TSMonitorSPEig",howoften,&howoften,NULL);CHKERRQ(ierr); 196 ierr = TSMonitorSPEigCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&ctx);CHKERRQ(ierr); 197 ierr = TSMonitorSet(ts,TSMonitorSPEig,ctx,(PetscErrorCode (*)(void**))TSMonitorSPEigCtxDestroy);CHKERRQ(ierr); 198 } 199 opt = PETSC_FALSE; 200 ierr = PetscOptionsName("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",&opt);CHKERRQ(ierr); 201 if (opt) { 202 TSMonitorDrawCtx ctx; 203 PetscInt howoften = 1; 204 205 ierr = PetscOptionsInt("-ts_monitor_draw_solution","Monitor solution graphically","TSMonitorDrawSolution",howoften,&howoften,NULL);CHKERRQ(ierr); 206 ierr = TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&ctx);CHKERRQ(ierr); 207 ierr = TSMonitorSet(ts,TSMonitorDrawSolution,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);CHKERRQ(ierr); 208 } 209 opt = PETSC_FALSE; 210 ierr = PetscOptionsName("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",&opt);CHKERRQ(ierr); 211 if (opt) { 212 TSMonitorDrawCtx ctx; 213 PetscReal bounds[4]; 214 PetscInt n = 4; 215 PetscDraw draw; 216 217 ierr = PetscOptionsRealArray("-ts_monitor_draw_solution_phase","Monitor solution graphically","TSMonitorDrawSolutionPhase",bounds,&n,NULL);CHKERRQ(ierr); 218 if (n != 4) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Must provide bounding box of phase field"); 219 ierr = TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,1,&ctx);CHKERRQ(ierr); 220 ierr = PetscViewerDrawGetDraw(ctx->viewer,0,&draw);CHKERRQ(ierr); 221 ierr = PetscDrawClear(draw);CHKERRQ(ierr); 222 ierr = PetscDrawAxisCreate(draw,&ctx->axis);CHKERRQ(ierr); 223 ierr = PetscDrawAxisSetLimits(ctx->axis,bounds[0],bounds[2],bounds[1],bounds[3]);CHKERRQ(ierr); 224 ierr = PetscDrawAxisSetLabels(ctx->axis,"Phase Diagram","Variable 1","Variable 2");CHKERRQ(ierr); 225 ierr = PetscDrawAxisDraw(ctx->axis);CHKERRQ(ierr); 226 /* ierr = PetscDrawSetCoordinates(draw,bounds[0],bounds[1],bounds[2],bounds[3]);CHKERRQ(ierr); */ 227 ierr = TSMonitorSet(ts,TSMonitorDrawSolutionPhase,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);CHKERRQ(ierr); 228 } 229 opt = PETSC_FALSE; 230 ierr = PetscOptionsName("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",&opt);CHKERRQ(ierr); 231 if (opt) { 232 TSMonitorDrawCtx ctx; 233 PetscInt howoften = 1; 234 235 ierr = PetscOptionsInt("-ts_monitor_draw_error","Monitor error graphically","TSMonitorDrawError",howoften,&howoften,NULL);CHKERRQ(ierr); 236 ierr = TSMonitorDrawCtxCreate(PetscObjectComm((PetscObject)ts),0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&ctx);CHKERRQ(ierr); 237 ierr = TSMonitorSet(ts,TSMonitorDrawError,ctx,(PetscErrorCode (*)(void**))TSMonitorDrawCtxDestroy);CHKERRQ(ierr); 238 } 239 opt = PETSC_FALSE; 240 ierr = PetscOptionsString("-ts_monitor_solution_binary","Save each solution to a binary file","TSMonitorSolutionBinary",0,monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 241 if (flg) { 242 PetscViewer ctx; 243 if (monfilename[0]) { 244 ierr = PetscViewerBinaryOpen(PetscObjectComm((PetscObject)ts),monfilename,FILE_MODE_WRITE,&ctx);CHKERRQ(ierr); 245 ierr = TSMonitorSet(ts,TSMonitorSolutionBinary,ctx,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 246 } else { 247 ctx = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)ts)); 248 ierr = TSMonitorSet(ts,TSMonitorSolutionBinary,ctx,(PetscErrorCode (*)(void**))NULL);CHKERRQ(ierr); 249 } 250 } 251 opt = PETSC_FALSE; 252 ierr = 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);CHKERRQ(ierr); 253 if (flg) { 254 const char *ptr,*ptr2; 255 char *filetemplate; 256 if (!monfilename[0]) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts"); 257 /* Do some cursory validation of the input. */ 258 ierr = PetscStrstr(monfilename,"%",(char**)&ptr);CHKERRQ(ierr); 259 if (!ptr) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"-ts_monitor_solution_vtk requires a file template, e.g. filename-%%03D.vts"); 260 for (ptr++; ptr && *ptr; ptr++) { 261 ierr = PetscStrchr("DdiouxX",*ptr,(char**)&ptr2);CHKERRQ(ierr); 262 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"); 263 if (ptr2) break; 264 } 265 ierr = PetscStrallocpy(monfilename,&filetemplate);CHKERRQ(ierr); 266 ierr = TSMonitorSet(ts,TSMonitorSolutionVTK,filetemplate,(PetscErrorCode (*)(void**))TSMonitorSolutionVTKDestroy);CHKERRQ(ierr); 267 } 268 269 ierr = PetscOptionsString("-ts_monitor_dmda_ray","Display a ray of the solution","None","y=0",dir,16,&flg);CHKERRQ(ierr); 270 if (flg) { 271 TSMonitorDMDARayCtx *rayctx; 272 int ray = 0; 273 DMDADirection ddir; 274 DM da; 275 PetscMPIInt rank; 276 277 if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir); 278 if (dir[0] == 'x') ddir = DMDA_X; 279 else if (dir[0] == 'y') ddir = DMDA_Y; 280 else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Unknown ray %s",dir); 281 sscanf(dir+2,"%d",&ray); 282 283 ierr = PetscInfo2(((PetscObject)ts),"Displaying DMDA ray %c = %D\n",dir[0],ray);CHKERRQ(ierr); 284 ierr = PetscNew(&rayctx);CHKERRQ(ierr); 285 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 286 ierr = DMDAGetRay(da,ddir,ray,&rayctx->ray,&rayctx->scatter);CHKERRQ(ierr); 287 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)ts),&rank);CHKERRQ(ierr); 288 if (!rank) { 289 ierr = PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,600,300,&rayctx->viewer);CHKERRQ(ierr); 290 } 291 rayctx->lgctx = NULL; 292 ierr = TSMonitorSet(ts,TSMonitorDMDARay,rayctx,TSMonitorDMDARayDestroy);CHKERRQ(ierr); 293 } 294 ierr = PetscOptionsString("-ts_monitor_lg_dmda_ray","Display a ray of the solution","None","x=0",dir,16,&flg);CHKERRQ(ierr); 295 if (flg) { 296 TSMonitorDMDARayCtx *rayctx; 297 int ray = 0; 298 DMDADirection ddir; 299 DM da; 300 PetscInt howoften = 1; 301 302 if (dir[1] != '=') SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Malformed ray %s", dir); 303 if (dir[0] == 'x') ddir = DMDA_X; 304 else if (dir[0] == 'y') ddir = DMDA_Y; 305 else SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Unknown ray direction %s", dir); 306 sscanf(dir+2, "%d", &ray); 307 308 ierr = PetscInfo2(((PetscObject) ts),"Displaying LG DMDA ray %c = %D\n", dir[0], ray);CHKERRQ(ierr); 309 ierr = PetscNew(&rayctx);CHKERRQ(ierr); 310 ierr = TSGetDM(ts, &da);CHKERRQ(ierr); 311 ierr = DMDAGetRay(da, ddir, ray, &rayctx->ray, &rayctx->scatter);CHKERRQ(ierr); 312 ierr = TSMonitorLGCtxCreate(PETSC_COMM_SELF,0,0,PETSC_DECIDE,PETSC_DECIDE,600,400,howoften,&rayctx->lgctx);CHKERRQ(ierr); 313 ierr = TSMonitorSet(ts, TSMonitorLGDMDARay, rayctx, TSMonitorDMDARayDestroy);CHKERRQ(ierr); 314 } 315 316 ierr = PetscOptionsName("-ts_monitor_envelope","Monitor maximum and minimum value of each component of the solution","TSMonitorEnvelope",&opt);CHKERRQ(ierr); 317 if (opt) { 318 TSMonitorEnvelopeCtx ctx; 319 320 ierr = TSMonitorEnvelopeCtxCreate(ts,&ctx);CHKERRQ(ierr); 321 ierr = TSMonitorSet(ts,TSMonitorEnvelope,ctx,(PetscErrorCode (*)(void**))TSMonitorEnvelopeCtxDestroy);CHKERRQ(ierr); 322 } 323 324 /* 325 This code is all wrong. One is creating objects inside the TSSetFromOptions() so if run with the options gui 326 will bleed memory. Also one is using a PetscOptionsBegin() inside a PetscOptionsBegin() 327 */ 328 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 329 ierr = TSAdaptSetFromOptions(PetscOptionsObject,adapt);CHKERRQ(ierr); 330 331 /* Handle specific TS options */ 332 if (ts->ops->setfromoptions) { 333 ierr = (*ts->ops->setfromoptions)(PetscOptionsObject,ts);CHKERRQ(ierr); 334 } 335 ierr = PetscOptionsEnd();CHKERRQ(ierr); 336 337 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 338 ierr = PetscObjectProcessOptionsHandlers((PetscObject)ts);CHKERRQ(ierr); 339 340 if (ts->trajectory) { 341 ierr = TSTrajectorySetFromOptions(ts->trajectory);CHKERRQ(ierr); 342 } 343 344 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 345 if (snes) { 346 if (ts->problem_type == TS_LINEAR) {ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr);} 347 ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 348 } 349 PetscFunctionReturn(0); 350 } 351 352 #undef __FUNCT__ 353 #define __FUNCT__ "TSSetSaveTrajectory" 354 /*@ 355 TSSetSaveTrajectory - Causes the TS to save its solutions as it iterates forward in time in a TSTrajectory object 356 357 Collective on TS 358 359 Input Parameters: 360 . ts - the TS context obtained from TSCreate() 361 362 363 Level: intermediate 364 365 .seealso: TSGetTrajectory(), TSAdjointSolve() 366 367 .keywords: TS, set, checkpoint, 368 @*/ 369 PetscErrorCode TSSetSaveTrajectory(TS ts) 370 { 371 PetscErrorCode ierr; 372 373 PetscFunctionBegin; 374 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 375 if (!ts->trajectory) { 376 ierr = TSTrajectoryCreate(PetscObjectComm((PetscObject)ts),&ts->trajectory);CHKERRQ(ierr); 377 ierr = TSTrajectorySetType(ts->trajectory,TSTRAJECTORYBASIC);CHKERRQ(ierr); 378 } 379 PetscFunctionReturn(0); 380 } 381 382 #undef __FUNCT__ 383 #define __FUNCT__ "TSComputeRHSJacobian" 384 /*@ 385 TSComputeRHSJacobian - Computes the Jacobian matrix that has been 386 set with TSSetRHSJacobian(). 387 388 Collective on TS and Vec 389 390 Input Parameters: 391 + ts - the TS context 392 . t - current timestep 393 - U - input vector 394 395 Output Parameters: 396 + A - Jacobian matrix 397 . B - optional preconditioning matrix 398 - flag - flag indicating matrix structure 399 400 Notes: 401 Most users should not need to explicitly call this routine, as it 402 is used internally within the nonlinear solvers. 403 404 See KSPSetOperators() for important information about setting the 405 flag parameter. 406 407 Level: developer 408 409 .keywords: SNES, compute, Jacobian, matrix 410 411 .seealso: TSSetRHSJacobian(), KSPSetOperators() 412 @*/ 413 PetscErrorCode TSComputeRHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B) 414 { 415 PetscErrorCode ierr; 416 PetscObjectState Ustate; 417 DM dm; 418 DMTS tsdm; 419 TSRHSJacobian rhsjacobianfunc; 420 void *ctx; 421 TSIJacobian ijacobianfunc; 422 423 PetscFunctionBegin; 424 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 425 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 426 PetscCheckSameComm(ts,1,U,3); 427 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 428 ierr = DMGetDMTS(dm,&tsdm);CHKERRQ(ierr); 429 ierr = DMTSGetRHSJacobian(dm,&rhsjacobianfunc,&ctx);CHKERRQ(ierr); 430 ierr = DMTSGetIJacobian(dm,&ijacobianfunc,NULL);CHKERRQ(ierr); 431 ierr = PetscObjectStateGet((PetscObject)U,&Ustate);CHKERRQ(ierr); 432 if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.X == U && ts->rhsjacobian.Xstate == Ustate))) { 433 PetscFunctionReturn(0); 434 } 435 436 if (!rhsjacobianfunc && !ijacobianfunc) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()"); 437 438 if (ts->rhsjacobian.reuse) { 439 ierr = MatShift(A,-ts->rhsjacobian.shift);CHKERRQ(ierr); 440 ierr = MatScale(A,1./ts->rhsjacobian.scale);CHKERRQ(ierr); 441 if (A != B) { 442 ierr = MatShift(B,-ts->rhsjacobian.shift);CHKERRQ(ierr); 443 ierr = MatScale(B,1./ts->rhsjacobian.scale);CHKERRQ(ierr); 444 } 445 ts->rhsjacobian.shift = 0; 446 ts->rhsjacobian.scale = 1.; 447 } 448 449 if (rhsjacobianfunc) { 450 ierr = PetscLogEventBegin(TS_JacobianEval,ts,U,A,B);CHKERRQ(ierr); 451 PetscStackPush("TS user Jacobian function"); 452 ierr = (*rhsjacobianfunc)(ts,t,U,A,B,ctx);CHKERRQ(ierr); 453 PetscStackPop; 454 ierr = PetscLogEventEnd(TS_JacobianEval,ts,U,A,B);CHKERRQ(ierr); 455 /* make sure user returned a correct Jacobian and preconditioner */ 456 PetscValidHeaderSpecific(A,MAT_CLASSID,4); 457 PetscValidHeaderSpecific(B,MAT_CLASSID,5); 458 } else { 459 ierr = MatZeroEntries(A);CHKERRQ(ierr); 460 if (A != B) {ierr = MatZeroEntries(B);CHKERRQ(ierr);} 461 } 462 ts->rhsjacobian.time = t; 463 ts->rhsjacobian.X = U; 464 ierr = PetscObjectStateGet((PetscObject)U,&ts->rhsjacobian.Xstate);CHKERRQ(ierr); 465 PetscFunctionReturn(0); 466 } 467 468 #undef __FUNCT__ 469 #define __FUNCT__ "TSComputeRHSFunction" 470 /*@ 471 TSComputeRHSFunction - Evaluates the right-hand-side function. 472 473 Collective on TS and Vec 474 475 Input Parameters: 476 + ts - the TS context 477 . t - current time 478 - U - state vector 479 480 Output Parameter: 481 . y - right hand side 482 483 Note: 484 Most users should not need to explicitly call this routine, as it 485 is used internally within the nonlinear solvers. 486 487 Level: developer 488 489 .keywords: TS, compute 490 491 .seealso: TSSetRHSFunction(), TSComputeIFunction() 492 @*/ 493 PetscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec U,Vec y) 494 { 495 PetscErrorCode ierr; 496 TSRHSFunction rhsfunction; 497 TSIFunction ifunction; 498 void *ctx; 499 DM dm; 500 501 PetscFunctionBegin; 502 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 503 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 504 PetscValidHeaderSpecific(y,VEC_CLASSID,4); 505 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 506 ierr = DMTSGetRHSFunction(dm,&rhsfunction,&ctx);CHKERRQ(ierr); 507 ierr = DMTSGetIFunction(dm,&ifunction,NULL);CHKERRQ(ierr); 508 509 if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()"); 510 511 ierr = PetscLogEventBegin(TS_FunctionEval,ts,U,y,0);CHKERRQ(ierr); 512 if (rhsfunction) { 513 PetscStackPush("TS user right-hand-side function"); 514 ierr = (*rhsfunction)(ts,t,U,y,ctx);CHKERRQ(ierr); 515 PetscStackPop; 516 } else { 517 ierr = VecZeroEntries(y);CHKERRQ(ierr); 518 } 519 520 ierr = PetscLogEventEnd(TS_FunctionEval,ts,U,y,0);CHKERRQ(ierr); 521 PetscFunctionReturn(0); 522 } 523 524 #undef __FUNCT__ 525 #define __FUNCT__ "TSComputeSolutionFunction" 526 /*@ 527 TSComputeSolutionFunction - Evaluates the solution function. 528 529 Collective on TS and Vec 530 531 Input Parameters: 532 + ts - the TS context 533 - t - current time 534 535 Output Parameter: 536 . U - the solution 537 538 Note: 539 Most users should not need to explicitly call this routine, as it 540 is used internally within the nonlinear solvers. 541 542 Level: developer 543 544 .keywords: TS, compute 545 546 .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction() 547 @*/ 548 PetscErrorCode TSComputeSolutionFunction(TS ts,PetscReal t,Vec U) 549 { 550 PetscErrorCode ierr; 551 TSSolutionFunction solutionfunction; 552 void *ctx; 553 DM dm; 554 555 PetscFunctionBegin; 556 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 557 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 558 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 559 ierr = DMTSGetSolutionFunction(dm,&solutionfunction,&ctx);CHKERRQ(ierr); 560 561 if (solutionfunction) { 562 PetscStackPush("TS user solution function"); 563 ierr = (*solutionfunction)(ts,t,U,ctx);CHKERRQ(ierr); 564 PetscStackPop; 565 } 566 PetscFunctionReturn(0); 567 } 568 #undef __FUNCT__ 569 #define __FUNCT__ "TSComputeForcingFunction" 570 /*@ 571 TSComputeForcingFunction - Evaluates the forcing function. 572 573 Collective on TS and Vec 574 575 Input Parameters: 576 + ts - the TS context 577 - t - current time 578 579 Output Parameter: 580 . U - the function value 581 582 Note: 583 Most users should not need to explicitly call this routine, as it 584 is used internally within the nonlinear solvers. 585 586 Level: developer 587 588 .keywords: TS, compute 589 590 .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction() 591 @*/ 592 PetscErrorCode TSComputeForcingFunction(TS ts,PetscReal t,Vec U) 593 { 594 PetscErrorCode ierr, (*forcing)(TS,PetscReal,Vec,void*); 595 void *ctx; 596 DM dm; 597 598 PetscFunctionBegin; 599 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 600 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 601 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 602 ierr = DMTSGetForcingFunction(dm,&forcing,&ctx);CHKERRQ(ierr); 603 604 if (forcing) { 605 PetscStackPush("TS user forcing function"); 606 ierr = (*forcing)(ts,t,U,ctx);CHKERRQ(ierr); 607 PetscStackPop; 608 } 609 PetscFunctionReturn(0); 610 } 611 612 #undef __FUNCT__ 613 #define __FUNCT__ "TSGetRHSVec_Private" 614 static PetscErrorCode TSGetRHSVec_Private(TS ts,Vec *Frhs) 615 { 616 Vec F; 617 PetscErrorCode ierr; 618 619 PetscFunctionBegin; 620 *Frhs = NULL; 621 ierr = TSGetIFunction(ts,&F,NULL,NULL);CHKERRQ(ierr); 622 if (!ts->Frhs) { 623 ierr = VecDuplicate(F,&ts->Frhs);CHKERRQ(ierr); 624 } 625 *Frhs = ts->Frhs; 626 PetscFunctionReturn(0); 627 } 628 629 #undef __FUNCT__ 630 #define __FUNCT__ "TSGetRHSMats_Private" 631 static PetscErrorCode TSGetRHSMats_Private(TS ts,Mat *Arhs,Mat *Brhs) 632 { 633 Mat A,B; 634 PetscErrorCode ierr; 635 636 PetscFunctionBegin; 637 if (Arhs) *Arhs = NULL; 638 if (Brhs) *Brhs = NULL; 639 ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr); 640 if (Arhs) { 641 if (!ts->Arhs) { 642 ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);CHKERRQ(ierr); 643 } 644 *Arhs = ts->Arhs; 645 } 646 if (Brhs) { 647 if (!ts->Brhs) { 648 if (A != B) { 649 ierr = MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ts->Brhs);CHKERRQ(ierr); 650 } else { 651 ts->Brhs = ts->Arhs; 652 ierr = PetscObjectReference((PetscObject)ts->Arhs);CHKERRQ(ierr); 653 } 654 } 655 *Brhs = ts->Brhs; 656 } 657 PetscFunctionReturn(0); 658 } 659 660 #undef __FUNCT__ 661 #define __FUNCT__ "TSComputeIFunction" 662 /*@ 663 TSComputeIFunction - Evaluates the DAE residual written in implicit form F(t,U,Udot)=0 664 665 Collective on TS and Vec 666 667 Input Parameters: 668 + ts - the TS context 669 . t - current time 670 . U - state vector 671 . Udot - time derivative of state vector 672 - imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate 673 674 Output Parameter: 675 . Y - right hand side 676 677 Note: 678 Most users should not need to explicitly call this routine, as it 679 is used internally within the nonlinear solvers. 680 681 If the user did did not write their equations in implicit form, this 682 function recasts them in implicit form. 683 684 Level: developer 685 686 .keywords: TS, compute 687 688 .seealso: TSSetIFunction(), TSComputeRHSFunction() 689 @*/ 690 PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec Y,PetscBool imex) 691 { 692 PetscErrorCode ierr; 693 TSIFunction ifunction; 694 TSRHSFunction rhsfunction; 695 void *ctx; 696 DM dm; 697 698 PetscFunctionBegin; 699 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 700 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 701 PetscValidHeaderSpecific(Udot,VEC_CLASSID,4); 702 PetscValidHeaderSpecific(Y,VEC_CLASSID,5); 703 704 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 705 ierr = DMTSGetIFunction(dm,&ifunction,&ctx);CHKERRQ(ierr); 706 ierr = DMTSGetRHSFunction(dm,&rhsfunction,NULL);CHKERRQ(ierr); 707 708 if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()"); 709 710 ierr = PetscLogEventBegin(TS_FunctionEval,ts,U,Udot,Y);CHKERRQ(ierr); 711 if (ifunction) { 712 PetscStackPush("TS user implicit function"); 713 ierr = (*ifunction)(ts,t,U,Udot,Y,ctx);CHKERRQ(ierr); 714 PetscStackPop; 715 } 716 if (imex) { 717 if (!ifunction) { 718 ierr = VecCopy(Udot,Y);CHKERRQ(ierr); 719 } 720 } else if (rhsfunction) { 721 if (ifunction) { 722 Vec Frhs; 723 ierr = TSGetRHSVec_Private(ts,&Frhs);CHKERRQ(ierr); 724 ierr = TSComputeRHSFunction(ts,t,U,Frhs);CHKERRQ(ierr); 725 ierr = VecAXPY(Y,-1,Frhs);CHKERRQ(ierr); 726 } else { 727 ierr = TSComputeRHSFunction(ts,t,U,Y);CHKERRQ(ierr); 728 ierr = VecAYPX(Y,-1,Udot);CHKERRQ(ierr); 729 } 730 } 731 ierr = PetscLogEventEnd(TS_FunctionEval,ts,U,Udot,Y);CHKERRQ(ierr); 732 PetscFunctionReturn(0); 733 } 734 735 #undef __FUNCT__ 736 #define __FUNCT__ "TSComputeIJacobian" 737 /*@ 738 TSComputeIJacobian - Evaluates the Jacobian of the DAE 739 740 Collective on TS and Vec 741 742 Input 743 Input Parameters: 744 + ts - the TS context 745 . t - current timestep 746 . U - state vector 747 . Udot - time derivative of state vector 748 . shift - shift to apply, see note below 749 - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate 750 751 Output Parameters: 752 + A - Jacobian matrix 753 . B - optional preconditioning matrix 754 - flag - flag indicating matrix structure 755 756 Notes: 757 If F(t,U,Udot)=0 is the DAE, the required Jacobian is 758 759 dF/dU + shift*dF/dUdot 760 761 Most users should not need to explicitly call this routine, as it 762 is used internally within the nonlinear solvers. 763 764 Level: developer 765 766 .keywords: TS, compute, Jacobian, matrix 767 768 .seealso: TSSetIJacobian() 769 @*/ 770 PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,PetscBool imex) 771 { 772 PetscErrorCode ierr; 773 TSIJacobian ijacobian; 774 TSRHSJacobian rhsjacobian; 775 DM dm; 776 void *ctx; 777 778 PetscFunctionBegin; 779 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 780 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 781 PetscValidHeaderSpecific(Udot,VEC_CLASSID,4); 782 PetscValidPointer(A,6); 783 PetscValidHeaderSpecific(A,MAT_CLASSID,6); 784 PetscValidPointer(B,7); 785 PetscValidHeaderSpecific(B,MAT_CLASSID,7); 786 787 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 788 ierr = DMTSGetIJacobian(dm,&ijacobian,&ctx);CHKERRQ(ierr); 789 ierr = DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);CHKERRQ(ierr); 790 791 if (!rhsjacobian && !ijacobian) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()"); 792 793 ierr = PetscLogEventBegin(TS_JacobianEval,ts,U,A,B);CHKERRQ(ierr); 794 if (ijacobian) { 795 PetscStackPush("TS user implicit Jacobian"); 796 ierr = (*ijacobian)(ts,t,U,Udot,shift,A,B,ctx);CHKERRQ(ierr); 797 PetscStackPop; 798 /* make sure user returned a correct Jacobian and preconditioner */ 799 PetscValidHeaderSpecific(A,MAT_CLASSID,4); 800 PetscValidHeaderSpecific(B,MAT_CLASSID,5); 801 } 802 if (imex) { 803 if (!ijacobian) { /* system was written as Udot = G(t,U) */ 804 ierr = MatZeroEntries(A);CHKERRQ(ierr); 805 ierr = MatShift(A,shift);CHKERRQ(ierr); 806 if (A != B) { 807 ierr = MatZeroEntries(B);CHKERRQ(ierr); 808 ierr = MatShift(B,shift);CHKERRQ(ierr); 809 } 810 } 811 } else { 812 Mat Arhs = NULL,Brhs = NULL; 813 if (rhsjacobian) { 814 if (ijacobian) { 815 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 816 } else { 817 ierr = TSGetIJacobian(ts,&Arhs,&Brhs,NULL,NULL);CHKERRQ(ierr); 818 } 819 ierr = TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);CHKERRQ(ierr); 820 } 821 if (Arhs == A) { /* No IJacobian, so we only have the RHS matrix */ 822 ts->rhsjacobian.scale = -1; 823 ts->rhsjacobian.shift = shift; 824 ierr = MatScale(A,-1);CHKERRQ(ierr); 825 ierr = MatShift(A,shift);CHKERRQ(ierr); 826 if (A != B) { 827 ierr = MatScale(B,-1);CHKERRQ(ierr); 828 ierr = MatShift(B,shift);CHKERRQ(ierr); 829 } 830 } else if (Arhs) { /* Both IJacobian and RHSJacobian */ 831 MatStructure axpy = DIFFERENT_NONZERO_PATTERN; 832 if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */ 833 ierr = MatZeroEntries(A);CHKERRQ(ierr); 834 ierr = MatShift(A,shift);CHKERRQ(ierr); 835 if (A != B) { 836 ierr = MatZeroEntries(B);CHKERRQ(ierr); 837 ierr = MatShift(B,shift);CHKERRQ(ierr); 838 } 839 } 840 ierr = MatAXPY(A,-1,Arhs,axpy);CHKERRQ(ierr); 841 if (A != B) { 842 ierr = MatAXPY(B,-1,Brhs,axpy);CHKERRQ(ierr); 843 } 844 } 845 } 846 ierr = PetscLogEventEnd(TS_JacobianEval,ts,U,A,B);CHKERRQ(ierr); 847 PetscFunctionReturn(0); 848 } 849 850 #undef __FUNCT__ 851 #define __FUNCT__ "TSSetRHSFunction" 852 /*@C 853 TSSetRHSFunction - Sets the routine for evaluating the function, 854 where U_t = G(t,u). 855 856 Logically Collective on TS 857 858 Input Parameters: 859 + ts - the TS context obtained from TSCreate() 860 . r - vector to put the computed right hand side (or NULL to have it created) 861 . f - routine for evaluating the right-hand-side function 862 - ctx - [optional] user-defined context for private data for the 863 function evaluation routine (may be NULL) 864 865 Calling sequence of func: 866 $ func (TS ts,PetscReal t,Vec u,Vec F,void *ctx); 867 868 + t - current timestep 869 . u - input vector 870 . F - function vector 871 - ctx - [optional] user-defined function context 872 873 Level: beginner 874 875 Notes: You must call this function or TSSetIFunction() to define your ODE. You cannot use this function when solving a DAE. 876 877 .keywords: TS, timestep, set, right-hand-side, function 878 879 .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSSetIFunction() 880 @*/ 881 PetscErrorCode TSSetRHSFunction(TS ts,Vec r,PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),void *ctx) 882 { 883 PetscErrorCode ierr; 884 SNES snes; 885 Vec ralloc = NULL; 886 DM dm; 887 888 PetscFunctionBegin; 889 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 890 if (r) PetscValidHeaderSpecific(r,VEC_CLASSID,2); 891 892 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 893 ierr = DMTSSetRHSFunction(dm,f,ctx);CHKERRQ(ierr); 894 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 895 if (!r && !ts->dm && ts->vec_sol) { 896 ierr = VecDuplicate(ts->vec_sol,&ralloc);CHKERRQ(ierr); 897 r = ralloc; 898 } 899 ierr = SNESSetFunction(snes,r,SNESTSFormFunction,ts);CHKERRQ(ierr); 900 ierr = VecDestroy(&ralloc);CHKERRQ(ierr); 901 PetscFunctionReturn(0); 902 } 903 904 #undef __FUNCT__ 905 #define __FUNCT__ "TSSetSolutionFunction" 906 /*@C 907 TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE 908 909 Logically Collective on TS 910 911 Input Parameters: 912 + ts - the TS context obtained from TSCreate() 913 . f - routine for evaluating the solution 914 - ctx - [optional] user-defined context for private data for the 915 function evaluation routine (may be NULL) 916 917 Calling sequence of func: 918 $ func (TS ts,PetscReal t,Vec u,void *ctx); 919 920 + t - current timestep 921 . u - output vector 922 - ctx - [optional] user-defined function context 923 924 Notes: 925 This routine is used for testing accuracy of time integration schemes when you already know the solution. 926 If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to 927 create closed-form solutions with non-physical forcing terms. 928 929 For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history. 930 931 Level: beginner 932 933 .keywords: TS, timestep, set, right-hand-side, function 934 935 .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetForcingFunction() 936 @*/ 937 PetscErrorCode TSSetSolutionFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx) 938 { 939 PetscErrorCode ierr; 940 DM dm; 941 942 PetscFunctionBegin; 943 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 944 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 945 ierr = DMTSSetSolutionFunction(dm,f,ctx);CHKERRQ(ierr); 946 PetscFunctionReturn(0); 947 } 948 949 #undef __FUNCT__ 950 #define __FUNCT__ "TSSetForcingFunction" 951 /*@C 952 TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE 953 954 Logically Collective on TS 955 956 Input Parameters: 957 + ts - the TS context obtained from TSCreate() 958 . f - routine for evaluating the forcing function 959 - ctx - [optional] user-defined context for private data for the 960 function evaluation routine (may be NULL) 961 962 Calling sequence of func: 963 $ func (TS ts,PetscReal t,Vec u,void *ctx); 964 965 + t - current timestep 966 . u - output vector 967 - ctx - [optional] user-defined function context 968 969 Notes: 970 This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to 971 create closed-form solutions with a non-physical forcing term. 972 973 For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history. 974 975 Level: beginner 976 977 .keywords: TS, timestep, set, right-hand-side, function 978 979 .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetSolutionFunction() 980 @*/ 981 PetscErrorCode TSSetForcingFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx) 982 { 983 PetscErrorCode ierr; 984 DM dm; 985 986 PetscFunctionBegin; 987 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 988 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 989 ierr = DMTSSetForcingFunction(dm,f,ctx);CHKERRQ(ierr); 990 PetscFunctionReturn(0); 991 } 992 993 #undef __FUNCT__ 994 #define __FUNCT__ "TSSetRHSJacobian" 995 /*@C 996 TSSetRHSJacobian - Sets the function to compute the Jacobian of G, 997 where U_t = G(U,t), as well as the location to store the matrix. 998 999 Logically Collective on TS 1000 1001 Input Parameters: 1002 + ts - the TS context obtained from TSCreate() 1003 . Amat - (approximate) Jacobian matrix 1004 . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat) 1005 . f - the Jacobian evaluation routine 1006 - ctx - [optional] user-defined context for private data for the 1007 Jacobian evaluation routine (may be NULL) 1008 1009 Calling sequence of f: 1010 $ func (TS ts,PetscReal t,Vec u,Mat A,Mat B,void *ctx); 1011 1012 + t - current timestep 1013 . u - input vector 1014 . Amat - (approximate) Jacobian matrix 1015 . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat) 1016 - ctx - [optional] user-defined context for matrix evaluation routine 1017 1018 1019 Level: beginner 1020 1021 .keywords: TS, timestep, set, right-hand-side, Jacobian 1022 1023 .seealso: SNESComputeJacobianDefaultColor(), TSSetRHSFunction(), TSRHSJacobianSetReuse(), TSSetIJacobian() 1024 1025 @*/ 1026 PetscErrorCode TSSetRHSJacobian(TS ts,Mat Amat,Mat Pmat,TSRHSJacobian f,void *ctx) 1027 { 1028 PetscErrorCode ierr; 1029 SNES snes; 1030 DM dm; 1031 TSIJacobian ijacobian; 1032 1033 PetscFunctionBegin; 1034 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1035 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 1036 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 1037 if (Amat) PetscCheckSameComm(ts,1,Amat,2); 1038 if (Pmat) PetscCheckSameComm(ts,1,Pmat,3); 1039 1040 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1041 ierr = DMTSSetRHSJacobian(dm,f,ctx);CHKERRQ(ierr); 1042 if (f == TSComputeRHSJacobianConstant) { 1043 /* Handle this case automatically for the user; otherwise user should call themselves. */ 1044 ierr = TSRHSJacobianSetReuse(ts,PETSC_TRUE);CHKERRQ(ierr); 1045 } 1046 ierr = DMTSGetIJacobian(dm,&ijacobian,NULL);CHKERRQ(ierr); 1047 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1048 if (!ijacobian) { 1049 ierr = SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1050 } 1051 if (Amat) { 1052 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 1053 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 1054 1055 ts->Arhs = Amat; 1056 } 1057 if (Pmat) { 1058 ierr = PetscObjectReference((PetscObject)Pmat);CHKERRQ(ierr); 1059 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 1060 1061 ts->Brhs = Pmat; 1062 } 1063 PetscFunctionReturn(0); 1064 } 1065 1066 1067 #undef __FUNCT__ 1068 #define __FUNCT__ "TSSetIFunction" 1069 /*@C 1070 TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved. 1071 1072 Logically Collective on TS 1073 1074 Input Parameters: 1075 + ts - the TS context obtained from TSCreate() 1076 . r - vector to hold the residual (or NULL to have it created internally) 1077 . f - the function evaluation routine 1078 - ctx - user-defined context for private data for the function evaluation routine (may be NULL) 1079 1080 Calling sequence of f: 1081 $ f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx); 1082 1083 + t - time at step/stage being solved 1084 . u - state vector 1085 . u_t - time derivative of state vector 1086 . F - function vector 1087 - ctx - [optional] user-defined context for matrix evaluation routine 1088 1089 Important: 1090 The user MUST call either this routine or TSSetRHSFunction() to define the ODE. When solving DAEs you must use this function. 1091 1092 Level: beginner 1093 1094 .keywords: TS, timestep, set, DAE, Jacobian 1095 1096 .seealso: TSSetRHSJacobian(), TSSetRHSFunction(), TSSetIJacobian() 1097 @*/ 1098 PetscErrorCode TSSetIFunction(TS ts,Vec res,TSIFunction f,void *ctx) 1099 { 1100 PetscErrorCode ierr; 1101 SNES snes; 1102 Vec resalloc = NULL; 1103 DM dm; 1104 1105 PetscFunctionBegin; 1106 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1107 if (res) PetscValidHeaderSpecific(res,VEC_CLASSID,2); 1108 1109 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1110 ierr = DMTSSetIFunction(dm,f,ctx);CHKERRQ(ierr); 1111 1112 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1113 if (!res && !ts->dm && ts->vec_sol) { 1114 ierr = VecDuplicate(ts->vec_sol,&resalloc);CHKERRQ(ierr); 1115 res = resalloc; 1116 } 1117 ierr = SNESSetFunction(snes,res,SNESTSFormFunction,ts);CHKERRQ(ierr); 1118 ierr = VecDestroy(&resalloc);CHKERRQ(ierr); 1119 PetscFunctionReturn(0); 1120 } 1121 1122 #undef __FUNCT__ 1123 #define __FUNCT__ "TSGetIFunction" 1124 /*@C 1125 TSGetIFunction - Returns the vector where the implicit residual is stored and the function/contex to compute it. 1126 1127 Not Collective 1128 1129 Input Parameter: 1130 . ts - the TS context 1131 1132 Output Parameter: 1133 + r - vector to hold residual (or NULL) 1134 . func - the function to compute residual (or NULL) 1135 - ctx - the function context (or NULL) 1136 1137 Level: advanced 1138 1139 .keywords: TS, nonlinear, get, function 1140 1141 .seealso: TSSetIFunction(), SNESGetFunction() 1142 @*/ 1143 PetscErrorCode TSGetIFunction(TS ts,Vec *r,TSIFunction *func,void **ctx) 1144 { 1145 PetscErrorCode ierr; 1146 SNES snes; 1147 DM dm; 1148 1149 PetscFunctionBegin; 1150 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1151 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1152 ierr = SNESGetFunction(snes,r,NULL,NULL);CHKERRQ(ierr); 1153 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1154 ierr = DMTSGetIFunction(dm,func,ctx);CHKERRQ(ierr); 1155 PetscFunctionReturn(0); 1156 } 1157 1158 #undef __FUNCT__ 1159 #define __FUNCT__ "TSGetRHSFunction" 1160 /*@C 1161 TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it. 1162 1163 Not Collective 1164 1165 Input Parameter: 1166 . ts - the TS context 1167 1168 Output Parameter: 1169 + r - vector to hold computed right hand side (or NULL) 1170 . func - the function to compute right hand side (or NULL) 1171 - ctx - the function context (or NULL) 1172 1173 Level: advanced 1174 1175 .keywords: TS, nonlinear, get, function 1176 1177 .seealso: TSSetRHSFunction(), SNESGetFunction() 1178 @*/ 1179 PetscErrorCode TSGetRHSFunction(TS ts,Vec *r,TSRHSFunction *func,void **ctx) 1180 { 1181 PetscErrorCode ierr; 1182 SNES snes; 1183 DM dm; 1184 1185 PetscFunctionBegin; 1186 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1187 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1188 ierr = SNESGetFunction(snes,r,NULL,NULL);CHKERRQ(ierr); 1189 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1190 ierr = DMTSGetRHSFunction(dm,func,ctx);CHKERRQ(ierr); 1191 PetscFunctionReturn(0); 1192 } 1193 1194 #undef __FUNCT__ 1195 #define __FUNCT__ "TSSetIJacobian" 1196 /*@C 1197 TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function 1198 provided with TSSetIFunction(). 1199 1200 Logically Collective on TS 1201 1202 Input Parameters: 1203 + ts - the TS context obtained from TSCreate() 1204 . Amat - (approximate) Jacobian matrix 1205 . Pmat - matrix used to compute preconditioner (usually the same as Amat) 1206 . f - the Jacobian evaluation routine 1207 - ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL) 1208 1209 Calling sequence of f: 1210 $ f(TS ts,PetscReal t,Vec U,Vec U_t,PetscReal a,Mat Amat,Mat Pmat,void *ctx); 1211 1212 + t - time at step/stage being solved 1213 . U - state vector 1214 . U_t - time derivative of state vector 1215 . a - shift 1216 . Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t 1217 . Pmat - matrix used for constructing preconditioner, usually the same as Amat 1218 - ctx - [optional] user-defined context for matrix evaluation routine 1219 1220 Notes: 1221 The matrices Amat and Pmat are exactly the matrices that are used by SNES for the nonlinear solve. 1222 1223 The matrix dF/dU + a*dF/dU_t you provide turns out to be 1224 the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved. 1225 The time integrator internally approximates U_t by W+a*U where the positive "shift" 1226 a and vector W depend on the integration method, step size, and past states. For example with 1227 the backward Euler method a = 1/dt and W = -a*U(previous timestep) so 1228 W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt 1229 1230 Level: beginner 1231 1232 .keywords: TS, timestep, DAE, Jacobian 1233 1234 .seealso: TSSetIFunction(), TSSetRHSJacobian(), SNESComputeJacobianDefaultColor(), SNESComputeJacobianDefault(), TSSetRHSFunction() 1235 1236 @*/ 1237 PetscErrorCode TSSetIJacobian(TS ts,Mat Amat,Mat Pmat,TSIJacobian f,void *ctx) 1238 { 1239 PetscErrorCode ierr; 1240 SNES snes; 1241 DM dm; 1242 1243 PetscFunctionBegin; 1244 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1245 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 1246 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 1247 if (Amat) PetscCheckSameComm(ts,1,Amat,2); 1248 if (Pmat) PetscCheckSameComm(ts,1,Pmat,3); 1249 1250 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1251 ierr = DMTSSetIJacobian(dm,f,ctx);CHKERRQ(ierr); 1252 1253 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1254 ierr = SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1255 PetscFunctionReturn(0); 1256 } 1257 1258 #undef __FUNCT__ 1259 #define __FUNCT__ "TSRHSJacobianSetReuse" 1260 /*@ 1261 TSRHSJacobianSetReuse - restore RHS Jacobian before re-evaluating. Without this flag, TS will change the sign and 1262 shift the RHS Jacobian for a finite-time-step implicit solve, in which case the user function will need to recompute 1263 the entire Jacobian. The reuse flag must be set if the evaluation function will assume that the matrix entries have 1264 not been changed by the TS. 1265 1266 Logically Collective 1267 1268 Input Arguments: 1269 + ts - TS context obtained from TSCreate() 1270 - reuse - PETSC_TRUE if the RHS Jacobian 1271 1272 Level: intermediate 1273 1274 .seealso: TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 1275 @*/ 1276 PetscErrorCode TSRHSJacobianSetReuse(TS ts,PetscBool reuse) 1277 { 1278 PetscFunctionBegin; 1279 ts->rhsjacobian.reuse = reuse; 1280 PetscFunctionReturn(0); 1281 } 1282 1283 #undef __FUNCT__ 1284 #define __FUNCT__ "TSLoad" 1285 /*@C 1286 TSLoad - Loads a KSP that has been stored in binary with KSPView(). 1287 1288 Collective on PetscViewer 1289 1290 Input Parameters: 1291 + newdm - the newly loaded TS, this needs to have been created with TSCreate() or 1292 some related function before a call to TSLoad(). 1293 - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() 1294 1295 Level: intermediate 1296 1297 Notes: 1298 The type is determined by the data in the file, any type set into the TS before this call is ignored. 1299 1300 Notes for advanced users: 1301 Most users should not need to know the details of the binary storage 1302 format, since TSLoad() and TSView() completely hide these details. 1303 But for anyone who's interested, the standard binary matrix storage 1304 format is 1305 .vb 1306 has not yet been determined 1307 .ve 1308 1309 .seealso: PetscViewerBinaryOpen(), TSView(), MatLoad(), VecLoad() 1310 @*/ 1311 PetscErrorCode TSLoad(TS ts, PetscViewer viewer) 1312 { 1313 PetscErrorCode ierr; 1314 PetscBool isbinary; 1315 PetscInt classid; 1316 char type[256]; 1317 DMTS sdm; 1318 DM dm; 1319 1320 PetscFunctionBegin; 1321 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1322 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1323 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1324 if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()"); 1325 1326 ierr = PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr); 1327 if (classid != TS_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Not TS next in file"); 1328 ierr = PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);CHKERRQ(ierr); 1329 ierr = TSSetType(ts, type);CHKERRQ(ierr); 1330 if (ts->ops->load) { 1331 ierr = (*ts->ops->load)(ts,viewer);CHKERRQ(ierr); 1332 } 1333 ierr = DMCreate(PetscObjectComm((PetscObject)ts),&dm);CHKERRQ(ierr); 1334 ierr = DMLoad(dm,viewer);CHKERRQ(ierr); 1335 ierr = TSSetDM(ts,dm);CHKERRQ(ierr); 1336 ierr = DMCreateGlobalVector(ts->dm,&ts->vec_sol);CHKERRQ(ierr); 1337 ierr = VecLoad(ts->vec_sol,viewer);CHKERRQ(ierr); 1338 ierr = DMGetDMTS(ts->dm,&sdm);CHKERRQ(ierr); 1339 ierr = DMTSLoad(sdm,viewer);CHKERRQ(ierr); 1340 PetscFunctionReturn(0); 1341 } 1342 1343 #include <petscdraw.h> 1344 #if defined(PETSC_HAVE_SAWS) 1345 #include <petscviewersaws.h> 1346 #endif 1347 #undef __FUNCT__ 1348 #define __FUNCT__ "TSView" 1349 /*@C 1350 TSView - Prints the TS data structure. 1351 1352 Collective on TS 1353 1354 Input Parameters: 1355 + ts - the TS context obtained from TSCreate() 1356 - viewer - visualization context 1357 1358 Options Database Key: 1359 . -ts_view - calls TSView() at end of TSStep() 1360 1361 Notes: 1362 The available visualization contexts include 1363 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 1364 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 1365 output where only the first processor opens 1366 the file. All other processors send their 1367 data to the first processor to print. 1368 1369 The user can open an alternative visualization context with 1370 PetscViewerASCIIOpen() - output to a specified file. 1371 1372 Level: beginner 1373 1374 .keywords: TS, timestep, view 1375 1376 .seealso: PetscViewerASCIIOpen() 1377 @*/ 1378 PetscErrorCode TSView(TS ts,PetscViewer viewer) 1379 { 1380 PetscErrorCode ierr; 1381 TSType type; 1382 PetscBool iascii,isstring,isundials,isbinary,isdraw; 1383 DMTS sdm; 1384 #if defined(PETSC_HAVE_SAWS) 1385 PetscBool issaws; 1386 #endif 1387 1388 PetscFunctionBegin; 1389 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1390 if (!viewer) { 1391 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts),&viewer);CHKERRQ(ierr); 1392 } 1393 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1394 PetscCheckSameComm(ts,1,viewer,2); 1395 1396 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1397 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 1398 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1399 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1400 #if defined(PETSC_HAVE_SAWS) 1401 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);CHKERRQ(ierr); 1402 #endif 1403 if (iascii) { 1404 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer);CHKERRQ(ierr); 1405 ierr = PetscViewerASCIIPrintf(viewer," maximum steps=%D\n",ts->max_steps);CHKERRQ(ierr); 1406 ierr = PetscViewerASCIIPrintf(viewer," maximum time=%g\n",(double)ts->max_time);CHKERRQ(ierr); 1407 if (ts->problem_type == TS_NONLINEAR) { 1408 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solver iterations=%D\n",ts->snes_its);CHKERRQ(ierr); 1409 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solve failures=%D\n",ts->num_snes_failures);CHKERRQ(ierr); 1410 } 1411 ierr = PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",ts->ksp_its);CHKERRQ(ierr); 1412 ierr = PetscViewerASCIIPrintf(viewer," total number of rejected steps=%D\n",ts->reject);CHKERRQ(ierr); 1413 ierr = DMGetDMTS(ts->dm,&sdm);CHKERRQ(ierr); 1414 ierr = DMTSView(sdm,viewer);CHKERRQ(ierr); 1415 if (ts->ops->view) { 1416 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1417 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1418 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1419 } 1420 } else if (isstring) { 1421 ierr = TSGetType(ts,&type);CHKERRQ(ierr); 1422 ierr = PetscViewerStringSPrintf(viewer," %-7.7s",type);CHKERRQ(ierr); 1423 } else if (isbinary) { 1424 PetscInt classid = TS_FILE_CLASSID; 1425 MPI_Comm comm; 1426 PetscMPIInt rank; 1427 char type[256]; 1428 1429 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 1430 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1431 if (!rank) { 1432 ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1433 ierr = PetscStrncpy(type,((PetscObject)ts)->type_name,256);CHKERRQ(ierr); 1434 ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr); 1435 } 1436 if (ts->ops->view) { 1437 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1438 } 1439 ierr = DMView(ts->dm,viewer);CHKERRQ(ierr); 1440 ierr = VecView(ts->vec_sol,viewer);CHKERRQ(ierr); 1441 ierr = DMGetDMTS(ts->dm,&sdm);CHKERRQ(ierr); 1442 ierr = DMTSView(sdm,viewer);CHKERRQ(ierr); 1443 } else if (isdraw) { 1444 PetscDraw draw; 1445 char str[36]; 1446 PetscReal x,y,bottom,h; 1447 1448 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1449 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 1450 ierr = PetscStrcpy(str,"TS: ");CHKERRQ(ierr); 1451 ierr = PetscStrcat(str,((PetscObject)ts)->type_name);CHKERRQ(ierr); 1452 ierr = PetscDrawBoxedString(draw,x,y,PETSC_DRAW_BLACK,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 1453 bottom = y - h; 1454 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr); 1455 if (ts->ops->view) { 1456 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1457 } 1458 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 1459 #if defined(PETSC_HAVE_SAWS) 1460 } else if (issaws) { 1461 PetscMPIInt rank; 1462 const char *name; 1463 1464 ierr = PetscObjectGetName((PetscObject)ts,&name);CHKERRQ(ierr); 1465 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 1466 if (!((PetscObject)ts)->amsmem && !rank) { 1467 char dir[1024]; 1468 1469 ierr = PetscObjectViewSAWs((PetscObject)ts,viewer);CHKERRQ(ierr); 1470 ierr = PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/time_step",name);CHKERRQ(ierr); 1471 PetscStackCallSAWs(SAWs_Register,(dir,&ts->steps,1,SAWs_READ,SAWs_INT)); 1472 ierr = PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/time",name);CHKERRQ(ierr); 1473 PetscStackCallSAWs(SAWs_Register,(dir,&ts->ptime,1,SAWs_READ,SAWs_DOUBLE)); 1474 } 1475 if (ts->ops->view) { 1476 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1477 } 1478 #endif 1479 } 1480 1481 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1482 ierr = PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);CHKERRQ(ierr); 1483 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1484 PetscFunctionReturn(0); 1485 } 1486 1487 1488 #undef __FUNCT__ 1489 #define __FUNCT__ "TSSetApplicationContext" 1490 /*@ 1491 TSSetApplicationContext - Sets an optional user-defined context for 1492 the timesteppers. 1493 1494 Logically Collective on TS 1495 1496 Input Parameters: 1497 + ts - the TS context obtained from TSCreate() 1498 - usrP - optional user context 1499 1500 Level: intermediate 1501 1502 .keywords: TS, timestep, set, application, context 1503 1504 .seealso: TSGetApplicationContext() 1505 @*/ 1506 PetscErrorCode TSSetApplicationContext(TS ts,void *usrP) 1507 { 1508 PetscFunctionBegin; 1509 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1510 ts->user = usrP; 1511 PetscFunctionReturn(0); 1512 } 1513 1514 #undef __FUNCT__ 1515 #define __FUNCT__ "TSGetApplicationContext" 1516 /*@ 1517 TSGetApplicationContext - Gets the user-defined context for the 1518 timestepper. 1519 1520 Not Collective 1521 1522 Input Parameter: 1523 . ts - the TS context obtained from TSCreate() 1524 1525 Output Parameter: 1526 . usrP - user context 1527 1528 Level: intermediate 1529 1530 .keywords: TS, timestep, get, application, context 1531 1532 .seealso: TSSetApplicationContext() 1533 @*/ 1534 PetscErrorCode TSGetApplicationContext(TS ts,void *usrP) 1535 { 1536 PetscFunctionBegin; 1537 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1538 *(void**)usrP = ts->user; 1539 PetscFunctionReturn(0); 1540 } 1541 1542 #undef __FUNCT__ 1543 #define __FUNCT__ "TSGetTimeStepNumber" 1544 /*@ 1545 TSGetTimeStepNumber - Gets the number of time steps completed. 1546 1547 Not Collective 1548 1549 Input Parameter: 1550 . ts - the TS context obtained from TSCreate() 1551 1552 Output Parameter: 1553 . iter - number of steps completed so far 1554 1555 Level: intermediate 1556 1557 .keywords: TS, timestep, get, iteration, number 1558 .seealso: TSGetTime(), TSGetTimeStep(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSSetPostStep() 1559 @*/ 1560 PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt *iter) 1561 { 1562 PetscFunctionBegin; 1563 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1564 PetscValidIntPointer(iter,2); 1565 *iter = ts->steps; 1566 PetscFunctionReturn(0); 1567 } 1568 1569 #undef __FUNCT__ 1570 #define __FUNCT__ "TSSetInitialTimeStep" 1571 /*@ 1572 TSSetInitialTimeStep - Sets the initial timestep to be used, 1573 as well as the initial time. 1574 1575 Logically Collective on TS 1576 1577 Input Parameters: 1578 + ts - the TS context obtained from TSCreate() 1579 . initial_time - the initial time 1580 - time_step - the size of the timestep 1581 1582 Level: intermediate 1583 1584 .seealso: TSSetTimeStep(), TSGetTimeStep() 1585 1586 .keywords: TS, set, initial, timestep 1587 @*/ 1588 PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step) 1589 { 1590 PetscErrorCode ierr; 1591 1592 PetscFunctionBegin; 1593 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1594 ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr); 1595 ierr = TSSetTime(ts,initial_time);CHKERRQ(ierr); 1596 PetscFunctionReturn(0); 1597 } 1598 1599 #undef __FUNCT__ 1600 #define __FUNCT__ "TSSetTimeStep" 1601 /*@ 1602 TSSetTimeStep - Allows one to reset the timestep at any time, 1603 useful for simple pseudo-timestepping codes. 1604 1605 Logically Collective on TS 1606 1607 Input Parameters: 1608 + ts - the TS context obtained from TSCreate() 1609 - time_step - the size of the timestep 1610 1611 Level: intermediate 1612 1613 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1614 1615 .keywords: TS, set, timestep 1616 @*/ 1617 PetscErrorCode TSSetTimeStep(TS ts,PetscReal time_step) 1618 { 1619 PetscFunctionBegin; 1620 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1621 PetscValidLogicalCollectiveReal(ts,time_step,2); 1622 ts->time_step = time_step; 1623 ts->time_step_orig = time_step; 1624 PetscFunctionReturn(0); 1625 } 1626 1627 #undef __FUNCT__ 1628 #define __FUNCT__ "TSSetExactFinalTime" 1629 /*@ 1630 TSSetExactFinalTime - Determines whether to adapt the final time step to 1631 match the exact final time, interpolate solution to the exact final time, 1632 or just return at the final time TS computed. 1633 1634 Logically Collective on TS 1635 1636 Input Parameter: 1637 + ts - the time-step context 1638 - eftopt - exact final time option 1639 1640 Level: beginner 1641 1642 .seealso: TSExactFinalTimeOption 1643 @*/ 1644 PetscErrorCode TSSetExactFinalTime(TS ts,TSExactFinalTimeOption eftopt) 1645 { 1646 PetscFunctionBegin; 1647 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1648 PetscValidLogicalCollectiveEnum(ts,eftopt,2); 1649 ts->exact_final_time = eftopt; 1650 PetscFunctionReturn(0); 1651 } 1652 1653 #undef __FUNCT__ 1654 #define __FUNCT__ "TSGetTimeStep" 1655 /*@ 1656 TSGetTimeStep - Gets the current timestep size. 1657 1658 Not Collective 1659 1660 Input Parameter: 1661 . ts - the TS context obtained from TSCreate() 1662 1663 Output Parameter: 1664 . dt - the current timestep size 1665 1666 Level: intermediate 1667 1668 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1669 1670 .keywords: TS, get, timestep 1671 @*/ 1672 PetscErrorCode TSGetTimeStep(TS ts,PetscReal *dt) 1673 { 1674 PetscFunctionBegin; 1675 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1676 PetscValidRealPointer(dt,2); 1677 *dt = ts->time_step; 1678 PetscFunctionReturn(0); 1679 } 1680 1681 #undef __FUNCT__ 1682 #define __FUNCT__ "TSGetSolution" 1683 /*@ 1684 TSGetSolution - Returns the solution at the present timestep. It 1685 is valid to call this routine inside the function that you are evaluating 1686 in order to move to the new timestep. This vector not changed until 1687 the solution at the next timestep has been calculated. 1688 1689 Not Collective, but Vec returned is parallel if TS is parallel 1690 1691 Input Parameter: 1692 . ts - the TS context obtained from TSCreate() 1693 1694 Output Parameter: 1695 . v - the vector containing the solution 1696 1697 Level: intermediate 1698 1699 .seealso: TSGetTimeStep() 1700 1701 .keywords: TS, timestep, get, solution 1702 @*/ 1703 PetscErrorCode TSGetSolution(TS ts,Vec *v) 1704 { 1705 PetscFunctionBegin; 1706 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1707 PetscValidPointer(v,2); 1708 *v = ts->vec_sol; 1709 PetscFunctionReturn(0); 1710 } 1711 1712 #undef __FUNCT__ 1713 #define __FUNCT__ "TSAdjointGetCostGradients" 1714 /*@ 1715 TSAdjointGetCostGradients - Returns the gradients from the TSAdjointSolve() 1716 1717 Not Collective, but Vec returned is parallel if TS is parallel 1718 1719 Input Parameter: 1720 . ts - the TS context obtained from TSCreate() 1721 1722 Output Parameter: 1723 + lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables 1724 - mu - vectors containing the gradients of the cost functions with respect to the problem parameters 1725 1726 Level: intermediate 1727 1728 .seealso: TSGetTimeStep() 1729 1730 .keywords: TS, timestep, get, sensitivity 1731 @*/ 1732 PetscErrorCode TSAdjointGetCostGradients(TS ts,PetscInt *numcost,Vec **lambda,Vec **mu) 1733 { 1734 PetscFunctionBegin; 1735 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1736 if (numcost) *numcost = ts->numcost; 1737 if (lambda) *lambda = ts->vecs_sensi; 1738 if (mu) *mu = ts->vecs_sensip; 1739 PetscFunctionReturn(0); 1740 } 1741 1742 /* ----- Routines to initialize and destroy a timestepper ---- */ 1743 #undef __FUNCT__ 1744 #define __FUNCT__ "TSSetProblemType" 1745 /*@ 1746 TSSetProblemType - Sets the type of problem to be solved. 1747 1748 Not collective 1749 1750 Input Parameters: 1751 + ts - The TS 1752 - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1753 .vb 1754 U_t - A U = 0 (linear) 1755 U_t - A(t) U = 0 (linear) 1756 F(t,U,U_t) = 0 (nonlinear) 1757 .ve 1758 1759 Level: beginner 1760 1761 .keywords: TS, problem type 1762 .seealso: TSSetUp(), TSProblemType, TS 1763 @*/ 1764 PetscErrorCode TSSetProblemType(TS ts, TSProblemType type) 1765 { 1766 PetscErrorCode ierr; 1767 1768 PetscFunctionBegin; 1769 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1770 ts->problem_type = type; 1771 if (type == TS_LINEAR) { 1772 SNES snes; 1773 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1774 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1775 } 1776 PetscFunctionReturn(0); 1777 } 1778 1779 #undef __FUNCT__ 1780 #define __FUNCT__ "TSGetProblemType" 1781 /*@C 1782 TSGetProblemType - Gets the type of problem to be solved. 1783 1784 Not collective 1785 1786 Input Parameter: 1787 . ts - The TS 1788 1789 Output Parameter: 1790 . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1791 .vb 1792 M U_t = A U 1793 M(t) U_t = A(t) U 1794 F(t,U,U_t) 1795 .ve 1796 1797 Level: beginner 1798 1799 .keywords: TS, problem type 1800 .seealso: TSSetUp(), TSProblemType, TS 1801 @*/ 1802 PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type) 1803 { 1804 PetscFunctionBegin; 1805 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1806 PetscValidIntPointer(type,2); 1807 *type = ts->problem_type; 1808 PetscFunctionReturn(0); 1809 } 1810 1811 #undef __FUNCT__ 1812 #define __FUNCT__ "TSSetUp" 1813 /*@ 1814 TSSetUp - Sets up the internal data structures for the later use 1815 of a timestepper. 1816 1817 Collective on TS 1818 1819 Input Parameter: 1820 . ts - the TS context obtained from TSCreate() 1821 1822 Notes: 1823 For basic use of the TS solvers the user need not explicitly call 1824 TSSetUp(), since these actions will automatically occur during 1825 the call to TSStep(). However, if one wishes to control this 1826 phase separately, TSSetUp() should be called after TSCreate() 1827 and optional routines of the form TSSetXXX(), but before TSStep(). 1828 1829 Level: advanced 1830 1831 .keywords: TS, timestep, setup 1832 1833 .seealso: TSCreate(), TSStep(), TSDestroy() 1834 @*/ 1835 PetscErrorCode TSSetUp(TS ts) 1836 { 1837 PetscErrorCode ierr; 1838 DM dm; 1839 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 1840 PetscErrorCode (*jac)(SNES,Vec,Mat,Mat,void*); 1841 TSIJacobian ijac; 1842 TSRHSJacobian rhsjac; 1843 1844 PetscFunctionBegin; 1845 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1846 if (ts->setupcalled) PetscFunctionReturn(0); 1847 1848 ts->total_steps = 0; 1849 if (!((PetscObject)ts)->type_name) { 1850 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 1851 } 1852 1853 if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first"); 1854 1855 1856 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 1857 1858 if (ts->rhsjacobian.reuse) { 1859 Mat Amat,Pmat; 1860 SNES snes; 1861 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1862 ierr = SNESGetJacobian(snes,&Amat,&Pmat,NULL,NULL);CHKERRQ(ierr); 1863 /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would 1864 * have displaced the RHS matrix */ 1865 if (Amat == ts->Arhs) { 1866 ierr = MatDuplicate(ts->Arhs,MAT_DO_NOT_COPY_VALUES,&Amat);CHKERRQ(ierr); 1867 ierr = SNESSetJacobian(snes,Amat,NULL,NULL,NULL);CHKERRQ(ierr); 1868 ierr = MatDestroy(&Amat);CHKERRQ(ierr); 1869 } 1870 if (Pmat == ts->Brhs) { 1871 ierr = MatDuplicate(ts->Brhs,MAT_DO_NOT_COPY_VALUES,&Pmat);CHKERRQ(ierr); 1872 ierr = SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);CHKERRQ(ierr); 1873 ierr = MatDestroy(&Pmat);CHKERRQ(ierr); 1874 } 1875 } 1876 if (ts->ops->setup) { 1877 ierr = (*ts->ops->setup)(ts);CHKERRQ(ierr); 1878 } 1879 1880 /* in the case where we've set a DMTSFunction or what have you, we need the default SNESFunction 1881 to be set right but can't do it elsewhere due to the overreliance on ctx=ts. 1882 */ 1883 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1884 ierr = DMSNESGetFunction(dm,&func,NULL);CHKERRQ(ierr); 1885 if (!func) { 1886 ierr =DMSNESSetFunction(dm,SNESTSFormFunction,ts);CHKERRQ(ierr); 1887 } 1888 /* if the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it. 1889 Otherwise, the SNES will use coloring internally to form the Jacobian. 1890 */ 1891 ierr = DMSNESGetJacobian(dm,&jac,NULL);CHKERRQ(ierr); 1892 ierr = DMTSGetIJacobian(dm,&ijac,NULL);CHKERRQ(ierr); 1893 ierr = DMTSGetRHSJacobian(dm,&rhsjac,NULL);CHKERRQ(ierr); 1894 if (!jac && (ijac || rhsjac)) { 1895 ierr = DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1896 } 1897 ts->setupcalled = PETSC_TRUE; 1898 PetscFunctionReturn(0); 1899 } 1900 1901 #undef __FUNCT__ 1902 #define __FUNCT__ "TSAdjointSetUp" 1903 /*@ 1904 TSAdjointSetUp - Sets up the internal data structures for the later use 1905 of an adjoint solver 1906 1907 Collective on TS 1908 1909 Input Parameter: 1910 . ts - the TS context obtained from TSCreate() 1911 1912 Notes: 1913 For basic use of the TS solvers the user need not explicitly call 1914 TSSetUp(), since these actions will automatically occur during 1915 the call to TSStep(). However, if one wishes to control this 1916 phase separately, TSSetUp() should be called after TSCreate() 1917 and optional routines of the form TSSetXXX(), but before TSStep(). 1918 1919 Level: advanced 1920 1921 .keywords: TS, timestep, setup 1922 1923 .seealso: TSCreate(), TSStep(), TSDestroy() 1924 @*/ 1925 PetscErrorCode TSAdjointSetUp(TS ts) 1926 { 1927 PetscErrorCode ierr; 1928 1929 PetscFunctionBegin; 1930 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1931 if (ts->adjointsetupcalled) PetscFunctionReturn(0); 1932 if (!ts->vecs_sensi) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSAdjointSetCostGradients() first"); 1933 1934 ierr = VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&ts->vecs_drdy);CHKERRQ(ierr); 1935 if (ts->vecs_sensip){ 1936 ierr = VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&ts->vecs_drdp);CHKERRQ(ierr); 1937 } 1938 1939 if (ts->ops->adjointsetup) { 1940 ierr = (*ts->ops->adjointsetup)(ts);CHKERRQ(ierr); 1941 } 1942 ts->adjointsetupcalled = PETSC_TRUE; 1943 PetscFunctionReturn(0); 1944 } 1945 1946 #undef __FUNCT__ 1947 #define __FUNCT__ "TSReset" 1948 /*@ 1949 TSReset - Resets a TS context and removes any allocated Vecs and Mats. 1950 1951 Collective on TS 1952 1953 Input Parameter: 1954 . ts - the TS context obtained from TSCreate() 1955 1956 Level: beginner 1957 1958 .keywords: TS, timestep, reset 1959 1960 .seealso: TSCreate(), TSSetup(), TSDestroy() 1961 @*/ 1962 PetscErrorCode TSReset(TS ts) 1963 { 1964 PetscErrorCode ierr; 1965 1966 PetscFunctionBegin; 1967 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1968 1969 if (ts->ops->reset) { 1970 ierr = (*ts->ops->reset)(ts);CHKERRQ(ierr); 1971 } 1972 if (ts->snes) {ierr = SNESReset(ts->snes);CHKERRQ(ierr);} 1973 1974 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 1975 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 1976 ierr = VecDestroy(&ts->Frhs);CHKERRQ(ierr); 1977 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1978 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 1979 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 1980 ierr = VecDestroyVecs(ts->nwork,&ts->work);CHKERRQ(ierr); 1981 ierr = VecDestroyVecs(ts->numcost,&ts->vecs_drdy);CHKERRQ(ierr); 1982 if (ts->vecs_drdp){ 1983 ierr = VecDestroyVecs(ts->numcost,&ts->vecs_drdp);CHKERRQ(ierr); 1984 } 1985 ts->vecs_sensi = NULL; 1986 ts->vecs_sensip = NULL; 1987 ierr = MatDestroy(&ts->Jacp);CHKERRQ(ierr); 1988 ierr = VecDestroy(&ts->vec_costintegral);CHKERRQ(ierr); 1989 ierr = VecDestroy(&ts->vec_costintegrand);CHKERRQ(ierr); 1990 ts->setupcalled = PETSC_FALSE; 1991 PetscFunctionReturn(0); 1992 } 1993 1994 #undef __FUNCT__ 1995 #define __FUNCT__ "TSDestroy" 1996 /*@ 1997 TSDestroy - Destroys the timestepper context that was created 1998 with TSCreate(). 1999 2000 Collective on TS 2001 2002 Input Parameter: 2003 . ts - the TS context obtained from TSCreate() 2004 2005 Level: beginner 2006 2007 .keywords: TS, timestepper, destroy 2008 2009 .seealso: TSCreate(), TSSetUp(), TSSolve() 2010 @*/ 2011 PetscErrorCode TSDestroy(TS *ts) 2012 { 2013 PetscErrorCode ierr; 2014 2015 PetscFunctionBegin; 2016 if (!*ts) PetscFunctionReturn(0); 2017 PetscValidHeaderSpecific((*ts),TS_CLASSID,1); 2018 if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; PetscFunctionReturn(0);} 2019 2020 ierr = TSReset((*ts));CHKERRQ(ierr); 2021 2022 /* if memory was published with SAWs then destroy it */ 2023 ierr = PetscObjectSAWsViewOff((PetscObject)*ts);CHKERRQ(ierr); 2024 if ((*ts)->ops->destroy) {ierr = (*(*ts)->ops->destroy)((*ts));CHKERRQ(ierr);} 2025 2026 ierr = TSTrajectoryDestroy(&(*ts)->trajectory);CHKERRQ(ierr); 2027 2028 ierr = TSAdaptDestroy(&(*ts)->adapt);CHKERRQ(ierr); 2029 if ((*ts)->event) { 2030 ierr = TSEventMonitorDestroy(&(*ts)->event);CHKERRQ(ierr); 2031 } 2032 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 2033 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 2034 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 2035 2036 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 2037 PetscFunctionReturn(0); 2038 } 2039 2040 #undef __FUNCT__ 2041 #define __FUNCT__ "TSGetSNES" 2042 /*@ 2043 TSGetSNES - Returns the SNES (nonlinear solver) associated with 2044 a TS (timestepper) context. Valid only for nonlinear problems. 2045 2046 Not Collective, but SNES is parallel if TS is parallel 2047 2048 Input Parameter: 2049 . ts - the TS context obtained from TSCreate() 2050 2051 Output Parameter: 2052 . snes - the nonlinear solver context 2053 2054 Notes: 2055 The user can then directly manipulate the SNES context to set various 2056 options, etc. Likewise, the user can then extract and manipulate the 2057 KSP, KSP, and PC contexts as well. 2058 2059 TSGetSNES() does not work for integrators that do not use SNES; in 2060 this case TSGetSNES() returns NULL in snes. 2061 2062 Level: beginner 2063 2064 .keywords: timestep, get, SNES 2065 @*/ 2066 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 2067 { 2068 PetscErrorCode ierr; 2069 2070 PetscFunctionBegin; 2071 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2072 PetscValidPointer(snes,2); 2073 if (!ts->snes) { 2074 ierr = SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);CHKERRQ(ierr); 2075 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 2076 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->snes);CHKERRQ(ierr); 2077 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 2078 if (ts->dm) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 2079 if (ts->problem_type == TS_LINEAR) { 2080 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 2081 } 2082 } 2083 *snes = ts->snes; 2084 PetscFunctionReturn(0); 2085 } 2086 2087 #undef __FUNCT__ 2088 #define __FUNCT__ "TSSetSNES" 2089 /*@ 2090 TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context 2091 2092 Collective 2093 2094 Input Parameter: 2095 + ts - the TS context obtained from TSCreate() 2096 - snes - the nonlinear solver context 2097 2098 Notes: 2099 Most users should have the TS created by calling TSGetSNES() 2100 2101 Level: developer 2102 2103 .keywords: timestep, set, SNES 2104 @*/ 2105 PetscErrorCode TSSetSNES(TS ts,SNES snes) 2106 { 2107 PetscErrorCode ierr; 2108 PetscErrorCode (*func)(SNES,Vec,Mat,Mat,void*); 2109 2110 PetscFunctionBegin; 2111 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2112 PetscValidHeaderSpecific(snes,SNES_CLASSID,2); 2113 ierr = PetscObjectReference((PetscObject)snes);CHKERRQ(ierr); 2114 ierr = SNESDestroy(&ts->snes);CHKERRQ(ierr); 2115 2116 ts->snes = snes; 2117 2118 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 2119 ierr = SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);CHKERRQ(ierr); 2120 if (func == SNESTSFormJacobian) { 2121 ierr = SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);CHKERRQ(ierr); 2122 } 2123 PetscFunctionReturn(0); 2124 } 2125 2126 #undef __FUNCT__ 2127 #define __FUNCT__ "TSGetKSP" 2128 /*@ 2129 TSGetKSP - Returns the KSP (linear solver) associated with 2130 a TS (timestepper) context. 2131 2132 Not Collective, but KSP is parallel if TS is parallel 2133 2134 Input Parameter: 2135 . ts - the TS context obtained from TSCreate() 2136 2137 Output Parameter: 2138 . ksp - the nonlinear solver context 2139 2140 Notes: 2141 The user can then directly manipulate the KSP context to set various 2142 options, etc. Likewise, the user can then extract and manipulate the 2143 KSP and PC contexts as well. 2144 2145 TSGetKSP() does not work for integrators that do not use KSP; 2146 in this case TSGetKSP() returns NULL in ksp. 2147 2148 Level: beginner 2149 2150 .keywords: timestep, get, KSP 2151 @*/ 2152 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 2153 { 2154 PetscErrorCode ierr; 2155 SNES snes; 2156 2157 PetscFunctionBegin; 2158 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2159 PetscValidPointer(ksp,2); 2160 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 2161 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 2162 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2163 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 2164 PetscFunctionReturn(0); 2165 } 2166 2167 /* ----------- Routines to set solver parameters ---------- */ 2168 2169 #undef __FUNCT__ 2170 #define __FUNCT__ "TSGetDuration" 2171 /*@ 2172 TSGetDuration - Gets the maximum number of timesteps to use and 2173 maximum time for iteration. 2174 2175 Not Collective 2176 2177 Input Parameters: 2178 + ts - the TS context obtained from TSCreate() 2179 . maxsteps - maximum number of iterations to use, or NULL 2180 - maxtime - final time to iterate to, or NULL 2181 2182 Level: intermediate 2183 2184 .keywords: TS, timestep, get, maximum, iterations, time 2185 @*/ 2186 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 2187 { 2188 PetscFunctionBegin; 2189 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2190 if (maxsteps) { 2191 PetscValidIntPointer(maxsteps,2); 2192 *maxsteps = ts->max_steps; 2193 } 2194 if (maxtime) { 2195 PetscValidScalarPointer(maxtime,3); 2196 *maxtime = ts->max_time; 2197 } 2198 PetscFunctionReturn(0); 2199 } 2200 2201 #undef __FUNCT__ 2202 #define __FUNCT__ "TSSetDuration" 2203 /*@ 2204 TSSetDuration - Sets the maximum number of timesteps to use and 2205 maximum time for iteration. 2206 2207 Logically Collective on TS 2208 2209 Input Parameters: 2210 + ts - the TS context obtained from TSCreate() 2211 . maxsteps - maximum number of iterations to use 2212 - maxtime - final time to iterate to 2213 2214 Options Database Keys: 2215 . -ts_max_steps <maxsteps> - Sets maxsteps 2216 . -ts_final_time <maxtime> - Sets maxtime 2217 2218 Notes: 2219 The default maximum number of iterations is 5000. Default time is 5.0 2220 2221 Level: intermediate 2222 2223 .keywords: TS, timestep, set, maximum, iterations 2224 2225 .seealso: TSSetExactFinalTime() 2226 @*/ 2227 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 2228 { 2229 PetscFunctionBegin; 2230 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2231 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 2232 PetscValidLogicalCollectiveReal(ts,maxtime,2); 2233 if (maxsteps >= 0) ts->max_steps = maxsteps; 2234 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 2235 PetscFunctionReturn(0); 2236 } 2237 2238 #undef __FUNCT__ 2239 #define __FUNCT__ "TSSetSolution" 2240 /*@ 2241 TSSetSolution - Sets the initial solution vector 2242 for use by the TS routines. 2243 2244 Logically Collective on TS and Vec 2245 2246 Input Parameters: 2247 + ts - the TS context obtained from TSCreate() 2248 - u - the solution vector 2249 2250 Level: beginner 2251 2252 .keywords: TS, timestep, set, solution, initial conditions 2253 @*/ 2254 PetscErrorCode TSSetSolution(TS ts,Vec u) 2255 { 2256 PetscErrorCode ierr; 2257 DM dm; 2258 2259 PetscFunctionBegin; 2260 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2261 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2262 ierr = PetscObjectReference((PetscObject)u);CHKERRQ(ierr); 2263 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 2264 2265 ts->vec_sol = u; 2266 2267 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2268 ierr = DMShellSetGlobalVector(dm,u);CHKERRQ(ierr); 2269 PetscFunctionReturn(0); 2270 } 2271 2272 #undef __FUNCT__ 2273 #define __FUNCT__ "TSAdjointSetSteps" 2274 /*@ 2275 TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time 2276 2277 Logically Collective on TS 2278 2279 Input Parameters: 2280 + ts - the TS context obtained from TSCreate() 2281 . steps - number of steps to use 2282 2283 Level: intermediate 2284 2285 Notes: Normally one does not call this and TSAdjointSolve() integrates back to the original timestep. One can call this 2286 so as to integrate back to less than the original timestep 2287 2288 .keywords: TS, timestep, set, maximum, iterations 2289 2290 .seealso: TSSetExactFinalTime() 2291 @*/ 2292 PetscErrorCode TSAdjointSetSteps(TS ts,PetscInt steps) 2293 { 2294 PetscFunctionBegin; 2295 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2296 PetscValidLogicalCollectiveInt(ts,steps,2); 2297 if (steps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back a negative number of steps"); 2298 if (steps > ts->total_steps) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back more than the total number of forward steps"); 2299 ts->adjoint_max_steps = steps; 2300 PetscFunctionReturn(0); 2301 } 2302 2303 #undef __FUNCT__ 2304 #define __FUNCT__ "TSAdjointSetCostGradients" 2305 /*@ 2306 TSAdjointSetCostGradients - Sets the initial value of the gradients of the cost function w.r.t. initial conditions and w.r.t. the problem parameters 2307 for use by the TSAdjoint routines. 2308 2309 Logically Collective on TS and Vec 2310 2311 Input Parameters: 2312 + ts - the TS context obtained from TSCreate() 2313 . lambda - gradients with respect to the initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector 2314 - mu - gradients with respect to the parameters, the number of entries in these vectors is the same as the number of parameters 2315 2316 Level: beginner 2317 2318 Notes: the entries in these vectors must be correctly initialized with the values lamda_i = df/dy|finaltime mu_i = df/dp|finaltime 2319 2320 .keywords: TS, timestep, set, sensitivity, initial conditions 2321 @*/ 2322 PetscErrorCode TSAdjointSetCostGradients(TS ts,PetscInt numcost,Vec *lambda,Vec *mu) 2323 { 2324 PetscFunctionBegin; 2325 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2326 PetscValidPointer(lambda,2); 2327 ts->vecs_sensi = lambda; 2328 ts->vecs_sensip = mu; 2329 ts->numcost = numcost; 2330 PetscFunctionReturn(0); 2331 } 2332 2333 #undef __FUNCT__ 2334 #define __FUNCT__ "TSAdjointSetRHSJacobian" 2335 /*@C 2336 TSAdjointSetRHSJacobian - Sets the function that computes the Jacobian of G w.r.t. the parameters p where y_t = G(y,p,t), as well as the location to store the matrix. 2337 2338 Logically Collective on TS 2339 2340 Input Parameters: 2341 + ts - The TS context obtained from TSCreate() 2342 - func - The function 2343 2344 Calling sequence of func: 2345 $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx); 2346 + t - current timestep 2347 . y - input vector (current ODE solution) 2348 . A - output matrix 2349 - ctx - [optional] user-defined function context 2350 2351 Level: intermediate 2352 2353 Notes: Amat has the same number of rows and the same row parallel layout as u, Amat has the same number of columns and parallel layout as p 2354 2355 .keywords: TS, sensitivity 2356 .seealso: 2357 @*/ 2358 PetscErrorCode TSAdjointSetRHSJacobian(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Mat,void*),void *ctx) 2359 { 2360 PetscErrorCode ierr; 2361 2362 PetscFunctionBegin; 2363 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2364 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 2365 2366 ts->rhsjacobianp = func; 2367 ts->rhsjacobianpctx = ctx; 2368 if(Amat) { 2369 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 2370 ierr = MatDestroy(&ts->Jacp);CHKERRQ(ierr); 2371 ts->Jacp = Amat; 2372 } 2373 PetscFunctionReturn(0); 2374 } 2375 2376 #undef __FUNCT__ 2377 #define __FUNCT__ "TSAdjointComputeRHSJacobian" 2378 /*@C 2379 TSAdjointComputeRHSJacobian - Runs the user-defined Jacobian function. 2380 2381 Collective on TS 2382 2383 Input Parameters: 2384 . ts - The TS context obtained from TSCreate() 2385 2386 Level: developer 2387 2388 .keywords: TS, sensitivity 2389 .seealso: TSAdjointSetRHSJacobian() 2390 @*/ 2391 PetscErrorCode TSAdjointComputeRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat) 2392 { 2393 PetscErrorCode ierr; 2394 2395 PetscFunctionBegin; 2396 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2397 PetscValidHeaderSpecific(X,VEC_CLASSID,3); 2398 PetscValidPointer(Amat,4); 2399 2400 PetscStackPush("TS user JacobianP function for sensitivity analysis"); 2401 ierr = (*ts->rhsjacobianp)(ts,t,X,Amat,ts->rhsjacobianpctx); CHKERRQ(ierr); 2402 PetscStackPop; 2403 PetscFunctionReturn(0); 2404 } 2405 2406 #undef __FUNCT__ 2407 #define __FUNCT__ "TSAdjointSetCostIntegrand" 2408 /*@C 2409 TSAdjointSetCostIntegrand - Sets the routine for evaluating the integral term in one or more cost functions 2410 2411 Logically Collective on TS 2412 2413 Input Parameters: 2414 + ts - the TS context obtained from TSCreate() 2415 . numcost - number of gradients to be computed, this is the number of cost functions 2416 . rf - routine for evaluating the integrand function 2417 . drdyf - function that computes the gradients of the r's with respect to y,NULL if not a function y 2418 . drdpf - function that computes the gradients of the r's with respect to p, NULL if not a function of p 2419 - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 2420 2421 Calling sequence of rf: 2422 $ rf(TS ts,PetscReal t,Vec y,Vec f[],void *ctx); 2423 2424 + t - current timestep 2425 . y - input vector 2426 . f - function result; one vector entry for each cost function 2427 - ctx - [optional] user-defined function context 2428 2429 Calling sequence of drdyf: 2430 $ PetscErroCode drdyf(TS ts,PetscReal t,Vec y,Vec *drdy,void *ctx); 2431 2432 Calling sequence of drdpf: 2433 $ PetscErroCode drdpf(TS ts,PetscReal t,Vec y,Vec *drdp,void *ctx); 2434 2435 Level: intermediate 2436 2437 Notes: For optimization there is generally a single cost function, numcost = 1. For sensitivities there may be multiple cost functions 2438 2439 .keywords: TS, sensitivity analysis, timestep, set, quadrature, function 2440 2441 .seealso: TSAdjointSetRHSJacobian(),TSAdjointGetCostGradients(), TSAdjointSetCostGradients() 2442 @*/ 2443 PetscErrorCode TSAdjointSetCostIntegrand(TS ts,PetscInt numcost, PetscErrorCode (*rf)(TS,PetscReal,Vec,Vec,void*), 2444 PetscErrorCode (*drdyf)(TS,PetscReal,Vec,Vec*,void*), 2445 PetscErrorCode (*drdpf)(TS,PetscReal,Vec,Vec*,void*),void *ctx) 2446 { 2447 PetscErrorCode ierr; 2448 2449 PetscFunctionBegin; 2450 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2451 if (!ts->numcost) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Call TSAdjointSetCostGradients() first so that the number of cost functions can be determined."); 2452 if (ts->numcost && ts->numcost!=numcost) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2rd parameter of TSAdjointSetCostIntegrand()) is inconsistent with the one set by TSAdjointSetCostGradients()"); 2453 2454 ierr = VecCreateSeq(PETSC_COMM_SELF,numcost,&ts->vec_costintegral);CHKERRQ(ierr); 2455 ierr = VecDuplicate(ts->vec_costintegral,&ts->vec_costintegrand);CHKERRQ(ierr); 2456 ts->costintegrand = rf; 2457 ts->costintegrandctx = ctx; 2458 ts->drdyfunction = drdyf; 2459 ts->drdpfunction = drdpf; 2460 PetscFunctionReturn(0); 2461 } 2462 2463 #undef __FUNCT__ 2464 #define __FUNCT__ "TSAdjointGetCostIntegral" 2465 /*@ 2466 TSAdjointGetCostIntegral - Returns the values of the integral term in the cost functions. 2467 It is valid to call the routine after a backward run. 2468 2469 Not Collective 2470 2471 Input Parameter: 2472 . ts - the TS context obtained from TSCreate() 2473 2474 Output Parameter: 2475 . v - the vector containing the integrals for each cost function 2476 2477 Level: intermediate 2478 2479 .seealso: TSAdjointSetCostIntegrand() 2480 2481 .keywords: TS, sensitivity analysis 2482 @*/ 2483 PetscErrorCode TSAdjointGetCostIntegral(TS ts,Vec *v) 2484 { 2485 PetscFunctionBegin; 2486 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2487 PetscValidPointer(v,2); 2488 *v = ts->vec_costintegral; 2489 PetscFunctionReturn(0); 2490 } 2491 2492 #undef __FUNCT__ 2493 #define __FUNCT__ "TSAdjointComputeCostIntegrand" 2494 /*@ 2495 TSAdjointComputeCostIntegrand - Evaluates the integral function in the cost functions. 2496 2497 Input Parameters: 2498 + ts - the TS context 2499 . t - current time 2500 - y - state vector, i.e. current solution 2501 2502 Output Parameter: 2503 . q - vector of size numcost to hold the outputs 2504 2505 Note: 2506 Most users should not need to explicitly call this routine, as it 2507 is used internally within the sensitivity analysis context. 2508 2509 Level: developer 2510 2511 .keywords: TS, compute 2512 2513 .seealso: TSAdjointSetCostIntegrand() 2514 @*/ 2515 PetscErrorCode TSAdjointComputeCostIntegrand(TS ts,PetscReal t,Vec y,Vec q) 2516 { 2517 PetscErrorCode ierr; 2518 2519 PetscFunctionBegin; 2520 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2521 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2522 PetscValidHeaderSpecific(q,VEC_CLASSID,4); 2523 2524 ierr = PetscLogEventBegin(TS_FunctionEval,ts,y,q,0);CHKERRQ(ierr); 2525 if (ts->costintegrand) { 2526 PetscStackPush("TS user integrand in the cost function"); 2527 ierr = (*ts->costintegrand)(ts,t,y,q,ts->costintegrandctx);CHKERRQ(ierr); 2528 PetscStackPop; 2529 } else { 2530 ierr = VecZeroEntries(q);CHKERRQ(ierr); 2531 } 2532 2533 ierr = PetscLogEventEnd(TS_FunctionEval,ts,y,q,0);CHKERRQ(ierr); 2534 PetscFunctionReturn(0); 2535 } 2536 2537 #undef __FUNCT__ 2538 #define __FUNCT__ "TSAdjointComputeDRDYFunction" 2539 /*@ 2540 TSAdjointComputeDRDYFunction - Runs the user-defined DRDY function. 2541 2542 Collective on TS 2543 2544 Input Parameters: 2545 . ts - The TS context obtained from TSCreate() 2546 2547 Notes: 2548 TSAdjointComputeDRDYFunction() is typically used for sensitivity implementation, 2549 so most users would not generally call this routine themselves. 2550 2551 Level: developer 2552 2553 .keywords: TS, sensitivity 2554 .seealso: TSAdjointComputeDRDYFunction() 2555 @*/ 2556 PetscErrorCode TSAdjointComputeDRDYFunction(TS ts,PetscReal t,Vec y,Vec *drdy) 2557 { 2558 PetscErrorCode ierr; 2559 2560 PetscFunctionBegin; 2561 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2562 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2563 2564 PetscStackPush("TS user DRDY function for sensitivity analysis"); 2565 ierr = (*ts->drdyfunction)(ts,t,y,drdy,ts->costintegrandctx); CHKERRQ(ierr); 2566 PetscStackPop; 2567 PetscFunctionReturn(0); 2568 } 2569 2570 #undef __FUNCT__ 2571 #define __FUNCT__ "TSAdjointComputeDRDPFunction" 2572 /*@ 2573 TSAdjointComputeDRDPFunction - Runs the user-defined DRDP function. 2574 2575 Collective on TS 2576 2577 Input Parameters: 2578 . ts - The TS context obtained from TSCreate() 2579 2580 Notes: 2581 TSDRDPFunction() is typically used for sensitivity implementation, 2582 so most users would not generally call this routine themselves. 2583 2584 Level: developer 2585 2586 .keywords: TS, sensitivity 2587 .seealso: TSAdjointSetDRDPFunction() 2588 @*/ 2589 PetscErrorCode TSAdjointComputeDRDPFunction(TS ts,PetscReal t,Vec y,Vec *drdp) 2590 { 2591 PetscErrorCode ierr; 2592 2593 PetscFunctionBegin; 2594 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2595 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2596 2597 PetscStackPush("TS user DRDP function for sensitivity analysis"); 2598 ierr = (*ts->drdpfunction)(ts,t,y,drdp,ts->costintegrandctx); CHKERRQ(ierr); 2599 PetscStackPop; 2600 PetscFunctionReturn(0); 2601 } 2602 2603 #undef __FUNCT__ 2604 #define __FUNCT__ "TSSetPreStep" 2605 /*@C 2606 TSSetPreStep - Sets the general-purpose function 2607 called once at the beginning of each time step. 2608 2609 Logically Collective on TS 2610 2611 Input Parameters: 2612 + ts - The TS context obtained from TSCreate() 2613 - func - The function 2614 2615 Calling sequence of func: 2616 . func (TS ts); 2617 2618 Level: intermediate 2619 2620 Note: 2621 If a step is rejected, TSStep() will call this routine again before each attempt. 2622 The last completed time step number can be queried using TSGetTimeStepNumber(), the 2623 size of the step being attempted can be obtained using TSGetTimeStep(). 2624 2625 .keywords: TS, timestep 2626 .seealso: TSSetPreStage(), TSSetPostStage(), TSSetPostStep(), TSStep() 2627 @*/ 2628 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 2629 { 2630 PetscFunctionBegin; 2631 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2632 ts->prestep = func; 2633 PetscFunctionReturn(0); 2634 } 2635 2636 #undef __FUNCT__ 2637 #define __FUNCT__ "TSPreStep" 2638 /*@ 2639 TSPreStep - Runs the user-defined pre-step function. 2640 2641 Collective on TS 2642 2643 Input Parameters: 2644 . ts - The TS context obtained from TSCreate() 2645 2646 Notes: 2647 TSPreStep() is typically used within time stepping implementations, 2648 so most users would not generally call this routine themselves. 2649 2650 Level: developer 2651 2652 .keywords: TS, timestep 2653 .seealso: TSSetPreStep(), TSPreStage(), TSPostStage(), TSPostStep() 2654 @*/ 2655 PetscErrorCode TSPreStep(TS ts) 2656 { 2657 PetscErrorCode ierr; 2658 2659 PetscFunctionBegin; 2660 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2661 if (ts->prestep) { 2662 PetscStackCallStandard((*ts->prestep),(ts)); 2663 } 2664 PetscFunctionReturn(0); 2665 } 2666 2667 #undef __FUNCT__ 2668 #define __FUNCT__ "TSSetPreStage" 2669 /*@C 2670 TSSetPreStage - Sets the general-purpose function 2671 called once at the beginning of each stage. 2672 2673 Logically Collective on TS 2674 2675 Input Parameters: 2676 + ts - The TS context obtained from TSCreate() 2677 - func - The function 2678 2679 Calling sequence of func: 2680 . PetscErrorCode func(TS ts, PetscReal stagetime); 2681 2682 Level: intermediate 2683 2684 Note: 2685 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2686 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2687 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2688 2689 .keywords: TS, timestep 2690 .seealso: TSSetPostStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2691 @*/ 2692 PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal)) 2693 { 2694 PetscFunctionBegin; 2695 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2696 ts->prestage = func; 2697 PetscFunctionReturn(0); 2698 } 2699 2700 #undef __FUNCT__ 2701 #define __FUNCT__ "TSSetPostStage" 2702 /*@C 2703 TSSetPostStage - Sets the general-purpose function 2704 called once at the end of each stage. 2705 2706 Logically Collective on TS 2707 2708 Input Parameters: 2709 + ts - The TS context obtained from TSCreate() 2710 - func - The function 2711 2712 Calling sequence of func: 2713 . PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y); 2714 2715 Level: intermediate 2716 2717 Note: 2718 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2719 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2720 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2721 2722 .keywords: TS, timestep 2723 .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2724 @*/ 2725 PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS,PetscReal,PetscInt,Vec*)) 2726 { 2727 PetscFunctionBegin; 2728 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2729 ts->poststage = func; 2730 PetscFunctionReturn(0); 2731 } 2732 2733 #undef __FUNCT__ 2734 #define __FUNCT__ "TSPreStage" 2735 /*@ 2736 TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage() 2737 2738 Collective on TS 2739 2740 Input Parameters: 2741 . ts - The TS context obtained from TSCreate() 2742 stagetime - The absolute time of the current stage 2743 2744 Notes: 2745 TSPreStage() is typically used within time stepping implementations, 2746 most users would not generally call this routine themselves. 2747 2748 Level: developer 2749 2750 .keywords: TS, timestep 2751 .seealso: TSPostStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 2752 @*/ 2753 PetscErrorCode TSPreStage(TS ts, PetscReal stagetime) 2754 { 2755 PetscErrorCode ierr; 2756 2757 PetscFunctionBegin; 2758 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2759 if (ts->prestage) { 2760 PetscStackCallStandard((*ts->prestage),(ts,stagetime)); 2761 } 2762 PetscFunctionReturn(0); 2763 } 2764 2765 #undef __FUNCT__ 2766 #define __FUNCT__ "TSPostStage" 2767 /*@ 2768 TSPostStage - Runs the user-defined post-stage function set using TSSetPostStage() 2769 2770 Collective on TS 2771 2772 Input Parameters: 2773 . ts - The TS context obtained from TSCreate() 2774 stagetime - The absolute time of the current stage 2775 stageindex - Stage number 2776 Y - Array of vectors (of size = total number 2777 of stages) with the stage solutions 2778 2779 Notes: 2780 TSPostStage() is typically used within time stepping implementations, 2781 most users would not generally call this routine themselves. 2782 2783 Level: developer 2784 2785 .keywords: TS, timestep 2786 .seealso: TSPreStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 2787 @*/ 2788 PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y) 2789 { 2790 PetscErrorCode ierr; 2791 2792 PetscFunctionBegin; 2793 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2794 if (ts->poststage) { 2795 PetscStackCallStandard((*ts->poststage),(ts,stagetime,stageindex,Y)); 2796 } 2797 PetscFunctionReturn(0); 2798 } 2799 2800 #undef __FUNCT__ 2801 #define __FUNCT__ "TSSetPostStep" 2802 /*@C 2803 TSSetPostStep - Sets the general-purpose function 2804 called once at the end of each time step. 2805 2806 Logically Collective on TS 2807 2808 Input Parameters: 2809 + ts - The TS context obtained from TSCreate() 2810 - func - The function 2811 2812 Calling sequence of func: 2813 $ func (TS ts); 2814 2815 Level: intermediate 2816 2817 .keywords: TS, timestep 2818 .seealso: TSSetPreStep(), TSSetPreStage(), TSGetTimeStep(), TSGetTimeStepNumber(), TSGetTime() 2819 @*/ 2820 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 2821 { 2822 PetscFunctionBegin; 2823 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2824 ts->poststep = func; 2825 PetscFunctionReturn(0); 2826 } 2827 2828 #undef __FUNCT__ 2829 #define __FUNCT__ "TSPostStep" 2830 /*@ 2831 TSPostStep - Runs the user-defined post-step function. 2832 2833 Collective on TS 2834 2835 Input Parameters: 2836 . ts - The TS context obtained from TSCreate() 2837 2838 Notes: 2839 TSPostStep() is typically used within time stepping implementations, 2840 so most users would not generally call this routine themselves. 2841 2842 Level: developer 2843 2844 .keywords: TS, timestep 2845 @*/ 2846 PetscErrorCode TSPostStep(TS ts) 2847 { 2848 PetscErrorCode ierr; 2849 2850 PetscFunctionBegin; 2851 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2852 if (ts->poststep) { 2853 PetscStackCallStandard((*ts->poststep),(ts)); 2854 } 2855 PetscFunctionReturn(0); 2856 } 2857 2858 /* ------------ Routines to set performance monitoring options ----------- */ 2859 2860 #undef __FUNCT__ 2861 #define __FUNCT__ "TSMonitorSet" 2862 /*@C 2863 TSMonitorSet - Sets an ADDITIONAL function that is to be used at every 2864 timestep to display the iteration's progress. 2865 2866 Logically Collective on TS 2867 2868 Input Parameters: 2869 + ts - the TS context obtained from TSCreate() 2870 . monitor - monitoring routine 2871 . mctx - [optional] user-defined context for private data for the 2872 monitor routine (use NULL if no context is desired) 2873 - monitordestroy - [optional] routine that frees monitor context 2874 (may be NULL) 2875 2876 Calling sequence of monitor: 2877 $ int monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx) 2878 2879 + ts - the TS context 2880 . steps - iteration number (after the final time step the monitor routine is called with a step of -1, this is at the final time which may have 2881 been interpolated to) 2882 . time - current time 2883 . u - current iterate 2884 - mctx - [optional] monitoring context 2885 2886 Notes: 2887 This routine adds an additional monitor to the list of monitors that 2888 already has been loaded. 2889 2890 Fortran notes: Only a single monitor function can be set for each TS object 2891 2892 Level: intermediate 2893 2894 .keywords: TS, timestep, set, monitor 2895 2896 .seealso: TSMonitorDefault(), TSMonitorCancel() 2897 @*/ 2898 PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**)) 2899 { 2900 PetscFunctionBegin; 2901 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2902 if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 2903 ts->monitor[ts->numbermonitors] = monitor; 2904 ts->monitordestroy[ts->numbermonitors] = mdestroy; 2905 ts->monitorcontext[ts->numbermonitors++] = (void*)mctx; 2906 PetscFunctionReturn(0); 2907 } 2908 2909 #undef __FUNCT__ 2910 #define __FUNCT__ "TSMonitorCancel" 2911 /*@C 2912 TSMonitorCancel - Clears all the monitors that have been set on a time-step object. 2913 2914 Logically Collective on TS 2915 2916 Input Parameters: 2917 . ts - the TS context obtained from TSCreate() 2918 2919 Notes: 2920 There is no way to remove a single, specific monitor. 2921 2922 Level: intermediate 2923 2924 .keywords: TS, timestep, set, monitor 2925 2926 .seealso: TSMonitorDefault(), TSMonitorSet() 2927 @*/ 2928 PetscErrorCode TSMonitorCancel(TS ts) 2929 { 2930 PetscErrorCode ierr; 2931 PetscInt i; 2932 2933 PetscFunctionBegin; 2934 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2935 for (i=0; i<ts->numbermonitors; i++) { 2936 if (ts->monitordestroy[i]) { 2937 ierr = (*ts->monitordestroy[i])(&ts->monitorcontext[i]);CHKERRQ(ierr); 2938 } 2939 } 2940 ts->numbermonitors = 0; 2941 PetscFunctionReturn(0); 2942 } 2943 2944 #undef __FUNCT__ 2945 #define __FUNCT__ "TSMonitorDefault" 2946 /*@ 2947 TSMonitorDefault - Sets the Default monitor 2948 2949 Level: intermediate 2950 2951 .keywords: TS, set, monitor 2952 2953 .seealso: TSMonitorDefault(), TSMonitorSet() 2954 @*/ 2955 PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy) 2956 { 2957 PetscErrorCode ierr; 2958 PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ts)); 2959 2960 PetscFunctionBegin; 2961 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2962 ierr = PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g\n",step,(double)ts->time_step,(double)ptime);CHKERRQ(ierr); 2963 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2964 PetscFunctionReturn(0); 2965 } 2966 2967 #undef __FUNCT__ 2968 #define __FUNCT__ "TSSetRetainStages" 2969 /*@ 2970 TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available. 2971 2972 Logically Collective on TS 2973 2974 Input Argument: 2975 . ts - time stepping context 2976 2977 Output Argument: 2978 . flg - PETSC_TRUE or PETSC_FALSE 2979 2980 Level: intermediate 2981 2982 .keywords: TS, set 2983 2984 .seealso: TSInterpolate(), TSSetPostStep() 2985 @*/ 2986 PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg) 2987 { 2988 PetscFunctionBegin; 2989 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2990 ts->retain_stages = flg; 2991 PetscFunctionReturn(0); 2992 } 2993 2994 #undef __FUNCT__ 2995 #define __FUNCT__ "TSInterpolate" 2996 /*@ 2997 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 2998 2999 Collective on TS 3000 3001 Input Argument: 3002 + ts - time stepping context 3003 - t - time to interpolate to 3004 3005 Output Argument: 3006 . U - state at given time 3007 3008 Notes: 3009 The user should call TSSetRetainStages() before taking a step in which interpolation will be requested. 3010 3011 Level: intermediate 3012 3013 Developer Notes: 3014 TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 3015 3016 .keywords: TS, set 3017 3018 .seealso: TSSetRetainStages(), TSSetPostStep() 3019 @*/ 3020 PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U) 3021 { 3022 PetscErrorCode ierr; 3023 3024 PetscFunctionBegin; 3025 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3026 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 3027 if (t < ts->ptime - ts->time_step_prev || t > ts->ptime) SETERRQ3(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Requested time %g not in last time steps [%g,%g]",t,(double)(ts->ptime-ts->time_step_prev),(double)ts->ptime); 3028 if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name); 3029 ierr = (*ts->ops->interpolate)(ts,t,U);CHKERRQ(ierr); 3030 PetscFunctionReturn(0); 3031 } 3032 3033 #undef __FUNCT__ 3034 #define __FUNCT__ "TSStep" 3035 /*@ 3036 TSStep - Steps one time step 3037 3038 Collective on TS 3039 3040 Input Parameter: 3041 . ts - the TS context obtained from TSCreate() 3042 3043 Level: developer 3044 3045 Notes: 3046 The public interface for the ODE/DAE solvers is TSSolve(), you should almost for sure be using that routine and not this routine. 3047 3048 The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may 3049 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 3050 3051 This may over-step the final time provided in TSSetDuration() depending on the time-step used. TSSolve() interpolates to exactly the 3052 time provided in TSSetDuration(). One can use TSInterpolate() to determine an interpolated solution within the final timestep. 3053 3054 .keywords: TS, timestep, solve 3055 3056 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSInterpolate() 3057 @*/ 3058 PetscErrorCode TSStep(TS ts) 3059 { 3060 DM dm; 3061 PetscErrorCode ierr; 3062 static PetscBool cite = PETSC_FALSE; 3063 3064 PetscFunctionBegin; 3065 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 3066 ierr = PetscCitationsRegister("@techreport{tspaper,\n" 3067 " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n" 3068 " author = {Shrirang Abhyankar and Jed Brown and Emil Constantinescu and Debojyoti Ghosh and Barry F. Smith},\n" 3069 " type = {Preprint},\n" 3070 " number = {ANL/MCS-P5061-0114},\n" 3071 " institution = {Argonne National Laboratory},\n" 3072 " year = {2014}\n}\n",&cite); 3073 3074 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 3075 ierr = TSSetUp(ts);CHKERRQ(ierr); 3076 3077 ts->reason = TS_CONVERGED_ITERATING; 3078 ts->ptime_prev = ts->ptime; 3079 ierr = DMSetOutputSequenceNumber(dm, ts->steps, ts->ptime);CHKERRQ(ierr); 3080 3081 if (!ts->ops->step) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSStep not implemented for type '%s'",((PetscObject)ts)->type_name); 3082 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 3083 ierr = (*ts->ops->step)(ts);CHKERRQ(ierr); 3084 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 3085 3086 ts->time_step_prev = ts->ptime - ts->ptime_prev; 3087 ierr = DMSetOutputSequenceNumber(dm, ts->steps, ts->ptime);CHKERRQ(ierr); 3088 3089 if (ts->reason < 0) { 3090 if (ts->errorifstepfailed) { 3091 if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s, increase -ts_max_snes_failures or make negative to attempt recovery",TSConvergedReasons[ts->reason]); 3092 else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 3093 } 3094 } else if (!ts->reason) { 3095 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 3096 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 3097 } 3098 ts->total_steps++; 3099 PetscFunctionReturn(0); 3100 } 3101 3102 #undef __FUNCT__ 3103 #define __FUNCT__ "TSAdjointStep" 3104 /*@ 3105 TSAdjointStep - Steps one time step 3106 3107 Collective on TS 3108 3109 Input Parameter: 3110 . ts - the TS context obtained from TSCreate() 3111 3112 Level: intermediate 3113 3114 Notes: 3115 The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may 3116 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 3117 3118 This may over-step the final time provided in TSSetDuration() depending on the time-step used. TSSolve() interpolates to exactly the 3119 time provided in TSSetDuration(). One can use TSInterpolate() to determine an interpolated solution within the final timestep. 3120 3121 .keywords: TS, timestep, solve 3122 3123 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSInterpolate() 3124 @*/ 3125 PetscErrorCode TSAdjointStep(TS ts) 3126 { 3127 DM dm; 3128 PetscErrorCode ierr; 3129 3130 PetscFunctionBegin; 3131 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 3132 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 3133 ierr = TSAdjointSetUp(ts);CHKERRQ(ierr); 3134 3135 ts->reason = TS_CONVERGED_ITERATING; 3136 ts->ptime_prev = ts->ptime; 3137 ierr = DMSetOutputSequenceNumber(dm, ts->steps, ts->ptime);CHKERRQ(ierr); 3138 ierr = VecViewFromOptions(ts->vec_sol, ((PetscObject) ts)->prefix, "-ts_view_solution");CHKERRQ(ierr); 3139 3140 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 3141 if (!ts->ops->adjointstep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed because the adjoint of %s has not been implemented, try other time stepping methods for adjoint sensitivity analysis",((PetscObject)ts)->type_name); 3142 ierr = (*ts->ops->adjointstep)(ts);CHKERRQ(ierr); 3143 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 3144 3145 ts->time_step_prev = ts->ptime - ts->ptime_prev; 3146 ierr = DMSetOutputSequenceNumber(dm, ts->steps, ts->ptime);CHKERRQ(ierr); 3147 3148 if (ts->reason < 0) { 3149 if (ts->errorifstepfailed) { 3150 if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) { 3151 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]); 3152 } else if (ts->reason == TS_DIVERGED_STEP_REJECTED) { 3153 SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s, increase -ts_max_reject or make negative to attempt recovery",TSConvergedReasons[ts->reason]); 3154 } else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 3155 } 3156 } else if (!ts->reason) { 3157 if (ts->steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS; 3158 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 3159 } 3160 ts->total_steps--; 3161 PetscFunctionReturn(0); 3162 } 3163 3164 #undef __FUNCT__ 3165 #define __FUNCT__ "TSEvaluateStep" 3166 /*@ 3167 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 3168 3169 Collective on TS 3170 3171 Input Arguments: 3172 + ts - time stepping context 3173 . order - desired order of accuracy 3174 - done - whether the step was evaluated at this order (pass NULL to generate an error if not available) 3175 3176 Output Arguments: 3177 . U - state at the end of the current step 3178 3179 Level: advanced 3180 3181 Notes: 3182 This function cannot be called until all stages have been evaluated. 3183 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. 3184 3185 .seealso: TSStep(), TSAdapt 3186 @*/ 3187 PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done) 3188 { 3189 PetscErrorCode ierr; 3190 3191 PetscFunctionBegin; 3192 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3193 PetscValidType(ts,1); 3194 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 3195 if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name); 3196 ierr = (*ts->ops->evaluatestep)(ts,order,U,done);CHKERRQ(ierr); 3197 PetscFunctionReturn(0); 3198 } 3199 3200 3201 #undef __FUNCT__ 3202 #define __FUNCT__ "TSSolve" 3203 /*@ 3204 TSSolve - Steps the requested number of timesteps. 3205 3206 Collective on TS 3207 3208 Input Parameter: 3209 + ts - the TS context obtained from TSCreate() 3210 - u - the solution vector (can be null if TSSetSolution() was used, otherwise must contain the initial conditions) 3211 3212 Level: beginner 3213 3214 Notes: 3215 The final time returned by this function may be different from the time of the internally 3216 held state accessible by TSGetSolution() and TSGetTime() because the method may have 3217 stepped over the final time. 3218 3219 .keywords: TS, timestep, solve 3220 3221 .seealso: TSCreate(), TSSetSolution(), TSStep() 3222 @*/ 3223 PetscErrorCode TSSolve(TS ts,Vec u) 3224 { 3225 Vec solution; 3226 PetscErrorCode ierr; 3227 3228 PetscFunctionBegin; 3229 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3230 if (u) PetscValidHeaderSpecific(u,VEC_CLASSID,2); 3231 if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE) { /* Need ts->vec_sol to be distinct so it is not overwritten when we interpolate at the end */ 3232 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 3233 if (!ts->vec_sol || u == ts->vec_sol) { 3234 ierr = VecDuplicate(u,&solution);CHKERRQ(ierr); 3235 ierr = TSSetSolution(ts,solution);CHKERRQ(ierr); 3236 ierr = VecDestroy(&solution);CHKERRQ(ierr); /* grant ownership */ 3237 } 3238 ierr = VecCopy(u,ts->vec_sol);CHKERRQ(ierr); 3239 } else if (u) { 3240 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 3241 } 3242 ierr = TSSetUp(ts);CHKERRQ(ierr); 3243 /* reset time step and iteration counters */ 3244 ts->steps = 0; 3245 ts->ksp_its = 0; 3246 ts->snes_its = 0; 3247 ts->num_snes_failures = 0; 3248 ts->reject = 0; 3249 ts->reason = TS_CONVERGED_ITERATING; 3250 3251 ierr = TSViewFromOptions(ts,NULL,"-ts_view_pre");CHKERRQ(ierr); 3252 3253 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 3254 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 3255 ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr); 3256 ts->solvetime = ts->ptime; 3257 } else { 3258 /* steps the requested number of timesteps. */ 3259 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 3260 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 3261 while (!ts->reason) { 3262 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 3263 ierr = TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 3264 ierr = TSStep(ts);CHKERRQ(ierr); 3265 if (ts->event) { 3266 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 3267 if (ts->event->status != TSEVENT_PROCESSING) { 3268 ierr = TSPostStep(ts);CHKERRQ(ierr); 3269 } 3270 } else { 3271 ierr = TSPostStep(ts);CHKERRQ(ierr); 3272 } 3273 } 3274 if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) { 3275 ierr = TSInterpolate(ts,ts->max_time,u);CHKERRQ(ierr); 3276 ts->solvetime = ts->max_time; 3277 solution = u; 3278 } else { 3279 if (u) {ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr);} 3280 ts->solvetime = ts->ptime; 3281 solution = ts->vec_sol; 3282 } 3283 ierr = TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->solvetime,solution);CHKERRQ(ierr); 3284 ierr = TSMonitor(ts,ts->steps,ts->solvetime,solution);CHKERRQ(ierr); 3285 ierr = VecViewFromOptions(solution, ((PetscObject) ts)->prefix, "-ts_view_solution");CHKERRQ(ierr); 3286 } 3287 3288 ierr = TSViewFromOptions(ts,NULL,"-ts_view");CHKERRQ(ierr); 3289 ierr = PetscObjectSAWsBlock((PetscObject)ts);CHKERRQ(ierr); 3290 if (ts->adjoint_solve) { 3291 ierr = TSAdjointSolve(ts);CHKERRQ(ierr); 3292 } 3293 PetscFunctionReturn(0); 3294 } 3295 3296 #undef __FUNCT__ 3297 #define __FUNCT__ "TSAdjointSolve" 3298 /*@ 3299 TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE 3300 3301 Collective on TS 3302 3303 Input Parameter: 3304 . ts - the TS context obtained from TSCreate() 3305 3306 Level: intermediate 3307 3308 Notes: 3309 This must be called after a call to TSSolve() that solves the forward problem 3310 3311 By default this will integrate back to the initial time, one can use TSAdjointSetSteps() to step back to a later time 3312 3313 .keywords: TS, timestep, solve 3314 3315 .seealso: TSCreate(), TSSetSolution(), TSStep() 3316 @*/ 3317 PetscErrorCode TSAdjointSolve(TS ts) 3318 { 3319 PetscErrorCode ierr; 3320 3321 PetscFunctionBegin; 3322 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3323 ierr = TSAdjointSetUp(ts);CHKERRQ(ierr); 3324 /* reset time step and iteration counters */ 3325 ts->steps = 0; 3326 ts->ksp_its = 0; 3327 ts->snes_its = 0; 3328 ts->num_snes_failures = 0; 3329 ts->reject = 0; 3330 ts->reason = TS_CONVERGED_ITERATING; 3331 3332 if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->total_steps; 3333 3334 if (ts->steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS; 3335 while (!ts->reason) { 3336 ierr = TSTrajectoryGet(ts->trajectory,ts,ts->adjoint_max_steps-ts->steps,ts->ptime);CHKERRQ(ierr); 3337 ierr = TSMonitor(ts,ts->adjoint_max_steps-ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 3338 ierr = TSAdjointStep(ts);CHKERRQ(ierr); 3339 if (ts->event) { 3340 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 3341 if (ts->event->status != TSEVENT_PROCESSING) { 3342 ierr = TSPostStep(ts);CHKERRQ(ierr); 3343 } 3344 } else { 3345 ierr = TSPostStep(ts);CHKERRQ(ierr); 3346 } 3347 } 3348 ts->solvetime = ts->ptime; 3349 PetscFunctionReturn(0); 3350 } 3351 3352 #undef __FUNCT__ 3353 #define __FUNCT__ "TSMonitor" 3354 /*@ 3355 TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet() 3356 3357 Collective on TS 3358 3359 Input Parameters: 3360 + ts - time stepping context obtained from TSCreate() 3361 . step - step number that has just completed 3362 . ptime - model time of the state 3363 - u - state at the current model time 3364 3365 Notes: 3366 TSMonitor() is typically used within the time stepping implementations. 3367 Users might call this function when using the TSStep() interface instead of TSSolve(). 3368 3369 Level: advanced 3370 3371 .keywords: TS, timestep 3372 @*/ 3373 PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u) 3374 { 3375 PetscErrorCode ierr; 3376 PetscInt i,n = ts->numbermonitors; 3377 3378 PetscFunctionBegin; 3379 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3380 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 3381 ierr = VecLockPush(u);CHKERRQ(ierr); 3382 for (i=0; i<n; i++) { 3383 ierr = (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);CHKERRQ(ierr); 3384 } 3385 ierr = VecLockPop(u);CHKERRQ(ierr); 3386 PetscFunctionReturn(0); 3387 } 3388 3389 /* ------------------------------------------------------------------------*/ 3390 #undef __FUNCT__ 3391 #define __FUNCT__ "TSMonitorLGCtxCreate" 3392 /*@C 3393 TSMonitorLGCtxCreate - Creates a line graph context for use with 3394 TS to monitor the solution process graphically in various ways 3395 3396 Collective on TS 3397 3398 Input Parameters: 3399 + host - the X display to open, or null for the local machine 3400 . label - the title to put in the title bar 3401 . x, y - the screen coordinates of the upper left coordinate of the window 3402 . m, n - the screen width and height in pixels 3403 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time 3404 3405 Output Parameter: 3406 . ctx - the context 3407 3408 Options Database Key: 3409 + -ts_monitor_lg_timestep - automatically sets line graph monitor 3410 . -ts_monitor_lg_solution - 3411 . -ts_monitor_lg_error - 3412 . -ts_monitor_lg_ksp_iterations - 3413 . -ts_monitor_lg_snes_iterations - 3414 - -lg_indicate_data_points <true,false> - indicate the data points (at each time step) on the plot; default is true 3415 3416 Notes: 3417 Use TSMonitorLGCtxDestroy() to destroy. 3418 3419 Level: intermediate 3420 3421 .keywords: TS, monitor, line graph, residual, seealso 3422 3423 .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 3424 3425 @*/ 3426 PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx) 3427 { 3428 PetscDraw win; 3429 PetscErrorCode ierr; 3430 3431 PetscFunctionBegin; 3432 ierr = PetscNew(ctx);CHKERRQ(ierr); 3433 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 3434 ierr = PetscDrawSetFromOptions(win);CHKERRQ(ierr); 3435 ierr = PetscDrawLGCreate(win,1,&(*ctx)->lg);CHKERRQ(ierr); 3436 ierr = PetscLogObjectParent((PetscObject)(*ctx)->lg,(PetscObject)win);CHKERRQ(ierr); 3437 ierr = PetscDrawLGIndicateDataPoints((*ctx)->lg,PETSC_TRUE);CHKERRQ(ierr); 3438 ierr = PetscDrawLGSetFromOptions((*ctx)->lg);CHKERRQ(ierr); 3439 (*ctx)->howoften = howoften; 3440 PetscFunctionReturn(0); 3441 } 3442 3443 #undef __FUNCT__ 3444 #define __FUNCT__ "TSMonitorLGTimeStep" 3445 PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx) 3446 { 3447 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 3448 PetscReal x = ptime,y; 3449 PetscErrorCode ierr; 3450 3451 PetscFunctionBegin; 3452 if (!step) { 3453 PetscDrawAxis axis; 3454 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 3455 ierr = PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time","Time step");CHKERRQ(ierr); 3456 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 3457 ierr = PetscDrawLGIndicateDataPoints(ctx->lg,PETSC_TRUE);CHKERRQ(ierr); 3458 } 3459 ierr = TSGetTimeStep(ts,&y);CHKERRQ(ierr); 3460 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 3461 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 3462 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 3463 } 3464 PetscFunctionReturn(0); 3465 } 3466 3467 #undef __FUNCT__ 3468 #define __FUNCT__ "TSMonitorLGCtxDestroy" 3469 /*@C 3470 TSMonitorLGCtxDestroy - Destroys a line graph context that was created 3471 with TSMonitorLGCtxCreate(). 3472 3473 Collective on TSMonitorLGCtx 3474 3475 Input Parameter: 3476 . ctx - the monitor context 3477 3478 Level: intermediate 3479 3480 .keywords: TS, monitor, line graph, destroy 3481 3482 .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep(); 3483 @*/ 3484 PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx) 3485 { 3486 PetscDraw draw; 3487 PetscErrorCode ierr; 3488 3489 PetscFunctionBegin; 3490 if ((*ctx)->transformdestroy) { 3491 ierr = ((*ctx)->transformdestroy)((*ctx)->transformctx);CHKERRQ(ierr); 3492 } 3493 ierr = PetscDrawLGGetDraw((*ctx)->lg,&draw);CHKERRQ(ierr); 3494 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 3495 ierr = PetscDrawLGDestroy(&(*ctx)->lg);CHKERRQ(ierr); 3496 ierr = PetscStrArrayDestroy(&(*ctx)->names);CHKERRQ(ierr); 3497 ierr = PetscStrArrayDestroy(&(*ctx)->displaynames);CHKERRQ(ierr); 3498 ierr = PetscFree((*ctx)->displayvariables);CHKERRQ(ierr); 3499 ierr = PetscFree((*ctx)->displayvalues);CHKERRQ(ierr); 3500 ierr = PetscFree(*ctx);CHKERRQ(ierr); 3501 PetscFunctionReturn(0); 3502 } 3503 3504 #undef __FUNCT__ 3505 #define __FUNCT__ "TSGetTime" 3506 /*@ 3507 TSGetTime - Gets the time of the most recently completed step. 3508 3509 Not Collective 3510 3511 Input Parameter: 3512 . ts - the TS context obtained from TSCreate() 3513 3514 Output Parameter: 3515 . t - the current time 3516 3517 Level: beginner 3518 3519 Note: 3520 When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(), 3521 TSSetPreStage(), TSSetPostStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated. 3522 3523 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 3524 3525 .keywords: TS, get, time 3526 @*/ 3527 PetscErrorCode TSGetTime(TS ts,PetscReal *t) 3528 { 3529 PetscFunctionBegin; 3530 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3531 PetscValidRealPointer(t,2); 3532 *t = ts->ptime; 3533 PetscFunctionReturn(0); 3534 } 3535 3536 #undef __FUNCT__ 3537 #define __FUNCT__ "TSGetPrevTime" 3538 /*@ 3539 TSGetPrevTime - Gets the starting time of the previously completed step. 3540 3541 Not Collective 3542 3543 Input Parameter: 3544 . ts - the TS context obtained from TSCreate() 3545 3546 Output Parameter: 3547 . t - the previous time 3548 3549 Level: beginner 3550 3551 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 3552 3553 .keywords: TS, get, time 3554 @*/ 3555 PetscErrorCode TSGetPrevTime(TS ts,PetscReal *t) 3556 { 3557 PetscFunctionBegin; 3558 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3559 PetscValidRealPointer(t,2); 3560 *t = ts->ptime_prev; 3561 PetscFunctionReturn(0); 3562 } 3563 3564 #undef __FUNCT__ 3565 #define __FUNCT__ "TSSetTime" 3566 /*@ 3567 TSSetTime - Allows one to reset the time. 3568 3569 Logically Collective on TS 3570 3571 Input Parameters: 3572 + ts - the TS context obtained from TSCreate() 3573 - time - the time 3574 3575 Level: intermediate 3576 3577 .seealso: TSGetTime(), TSSetDuration() 3578 3579 .keywords: TS, set, time 3580 @*/ 3581 PetscErrorCode TSSetTime(TS ts, PetscReal t) 3582 { 3583 PetscFunctionBegin; 3584 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3585 PetscValidLogicalCollectiveReal(ts,t,2); 3586 ts->ptime = t; 3587 PetscFunctionReturn(0); 3588 } 3589 3590 #undef __FUNCT__ 3591 #define __FUNCT__ "TSSetOptionsPrefix" 3592 /*@C 3593 TSSetOptionsPrefix - Sets the prefix used for searching for all 3594 TS options in the database. 3595 3596 Logically Collective on TS 3597 3598 Input Parameter: 3599 + ts - The TS context 3600 - prefix - The prefix to prepend to all option names 3601 3602 Notes: 3603 A hyphen (-) must NOT be given at the beginning of the prefix name. 3604 The first character of all runtime options is AUTOMATICALLY the 3605 hyphen. 3606 3607 Level: advanced 3608 3609 .keywords: TS, set, options, prefix, database 3610 3611 .seealso: TSSetFromOptions() 3612 3613 @*/ 3614 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 3615 { 3616 PetscErrorCode ierr; 3617 SNES snes; 3618 3619 PetscFunctionBegin; 3620 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3621 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3622 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3623 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 3624 PetscFunctionReturn(0); 3625 } 3626 3627 3628 #undef __FUNCT__ 3629 #define __FUNCT__ "TSAppendOptionsPrefix" 3630 /*@C 3631 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 3632 TS options in the database. 3633 3634 Logically Collective on TS 3635 3636 Input Parameter: 3637 + ts - The TS context 3638 - prefix - The prefix to prepend to all option names 3639 3640 Notes: 3641 A hyphen (-) must NOT be given at the beginning of the prefix name. 3642 The first character of all runtime options is AUTOMATICALLY the 3643 hyphen. 3644 3645 Level: advanced 3646 3647 .keywords: TS, append, options, prefix, database 3648 3649 .seealso: TSGetOptionsPrefix() 3650 3651 @*/ 3652 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 3653 { 3654 PetscErrorCode ierr; 3655 SNES snes; 3656 3657 PetscFunctionBegin; 3658 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3659 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3660 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3661 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 3662 PetscFunctionReturn(0); 3663 } 3664 3665 #undef __FUNCT__ 3666 #define __FUNCT__ "TSGetOptionsPrefix" 3667 /*@C 3668 TSGetOptionsPrefix - Sets the prefix used for searching for all 3669 TS options in the database. 3670 3671 Not Collective 3672 3673 Input Parameter: 3674 . ts - The TS context 3675 3676 Output Parameter: 3677 . prefix - A pointer to the prefix string used 3678 3679 Notes: On the fortran side, the user should pass in a string 'prifix' of 3680 sufficient length to hold the prefix. 3681 3682 Level: intermediate 3683 3684 .keywords: TS, get, options, prefix, database 3685 3686 .seealso: TSAppendOptionsPrefix() 3687 @*/ 3688 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 3689 { 3690 PetscErrorCode ierr; 3691 3692 PetscFunctionBegin; 3693 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3694 PetscValidPointer(prefix,2); 3695 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3696 PetscFunctionReturn(0); 3697 } 3698 3699 #undef __FUNCT__ 3700 #define __FUNCT__ "TSGetRHSJacobian" 3701 /*@C 3702 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 3703 3704 Not Collective, but parallel objects are returned if TS is parallel 3705 3706 Input Parameter: 3707 . ts - The TS context obtained from TSCreate() 3708 3709 Output Parameters: 3710 + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or NULL) 3711 . Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat (or NULL) 3712 . func - Function to compute the Jacobian of the RHS (or NULL) 3713 - ctx - User-defined context for Jacobian evaluation routine (or NULL) 3714 3715 Notes: You can pass in NULL for any return argument you do not need. 3716 3717 Level: intermediate 3718 3719 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3720 3721 .keywords: TS, timestep, get, matrix, Jacobian 3722 @*/ 3723 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx) 3724 { 3725 PetscErrorCode ierr; 3726 SNES snes; 3727 DM dm; 3728 3729 PetscFunctionBegin; 3730 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3731 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3732 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3733 ierr = DMTSGetRHSJacobian(dm,func,ctx);CHKERRQ(ierr); 3734 PetscFunctionReturn(0); 3735 } 3736 3737 #undef __FUNCT__ 3738 #define __FUNCT__ "TSGetIJacobian" 3739 /*@C 3740 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 3741 3742 Not Collective, but parallel objects are returned if TS is parallel 3743 3744 Input Parameter: 3745 . ts - The TS context obtained from TSCreate() 3746 3747 Output Parameters: 3748 + Amat - The (approximate) Jacobian of F(t,U,U_t) 3749 . Pmat - The matrix from which the preconditioner is constructed, often the same as Amat 3750 . f - The function to compute the matrices 3751 - ctx - User-defined context for Jacobian evaluation routine 3752 3753 Notes: You can pass in NULL for any return argument you do not need. 3754 3755 Level: advanced 3756 3757 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3758 3759 .keywords: TS, timestep, get, matrix, Jacobian 3760 @*/ 3761 PetscErrorCode TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx) 3762 { 3763 PetscErrorCode ierr; 3764 SNES snes; 3765 DM dm; 3766 3767 PetscFunctionBegin; 3768 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3769 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 3770 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3771 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3772 ierr = DMTSGetIJacobian(dm,f,ctx);CHKERRQ(ierr); 3773 PetscFunctionReturn(0); 3774 } 3775 3776 3777 #undef __FUNCT__ 3778 #define __FUNCT__ "TSMonitorDrawSolution" 3779 /*@C 3780 TSMonitorDrawSolution - Monitors progress of the TS solvers by calling 3781 VecView() for the solution at each timestep 3782 3783 Collective on TS 3784 3785 Input Parameters: 3786 + ts - the TS context 3787 . step - current time-step 3788 . ptime - current time 3789 - dummy - either a viewer or NULL 3790 3791 Options Database: 3792 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3793 3794 Notes: the initial solution and current solution are not display with a common axis scaling so generally the option -ts_monitor_draw_solution_initial 3795 will look bad 3796 3797 Level: intermediate 3798 3799 .keywords: TS, vector, monitor, view 3800 3801 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3802 @*/ 3803 PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3804 { 3805 PetscErrorCode ierr; 3806 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3807 PetscDraw draw; 3808 3809 PetscFunctionBegin; 3810 if (!step && ictx->showinitial) { 3811 if (!ictx->initialsolution) { 3812 ierr = VecDuplicate(u,&ictx->initialsolution);CHKERRQ(ierr); 3813 } 3814 ierr = VecCopy(u,ictx->initialsolution);CHKERRQ(ierr); 3815 } 3816 if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 3817 3818 if (ictx->showinitial) { 3819 PetscReal pause; 3820 ierr = PetscViewerDrawGetPause(ictx->viewer,&pause);CHKERRQ(ierr); 3821 ierr = PetscViewerDrawSetPause(ictx->viewer,0.0);CHKERRQ(ierr); 3822 ierr = VecView(ictx->initialsolution,ictx->viewer);CHKERRQ(ierr); 3823 ierr = PetscViewerDrawSetPause(ictx->viewer,pause);CHKERRQ(ierr); 3824 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);CHKERRQ(ierr); 3825 } 3826 ierr = VecView(u,ictx->viewer);CHKERRQ(ierr); 3827 if (ictx->showtimestepandtime) { 3828 PetscReal xl,yl,xr,yr,tw,w,h; 3829 char time[32]; 3830 size_t len; 3831 3832 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3833 ierr = PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);CHKERRQ(ierr); 3834 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3835 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3836 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3837 w = xl + .5*(xr - xl) - .5*len*tw; 3838 h = yl + .95*(yr - yl); 3839 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3840 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3841 } 3842 3843 if (ictx->showinitial) { 3844 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);CHKERRQ(ierr); 3845 } 3846 PetscFunctionReturn(0); 3847 } 3848 3849 #undef __FUNCT__ 3850 #define __FUNCT__ "TSMonitorDrawSolutionPhase" 3851 /*@C 3852 TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram 3853 3854 Collective on TS 3855 3856 Input Parameters: 3857 + ts - the TS context 3858 . step - current time-step 3859 . ptime - current time 3860 - dummy - either a viewer or NULL 3861 3862 Level: intermediate 3863 3864 .keywords: TS, vector, monitor, view 3865 3866 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3867 @*/ 3868 PetscErrorCode TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3869 { 3870 PetscErrorCode ierr; 3871 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3872 PetscDraw draw; 3873 MPI_Comm comm; 3874 PetscInt n; 3875 PetscMPIInt size; 3876 PetscReal xl,yl,xr,yr,tw,w,h; 3877 char time[32]; 3878 size_t len; 3879 const PetscScalar *U; 3880 3881 PetscFunctionBegin; 3882 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 3883 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3884 if (size != 1) SETERRQ(comm,PETSC_ERR_SUP,"Only allowed for sequential runs"); 3885 ierr = VecGetSize(u,&n);CHKERRQ(ierr); 3886 if (n != 2) SETERRQ(comm,PETSC_ERR_SUP,"Only for ODEs with two unknowns"); 3887 3888 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3889 3890 ierr = VecGetArrayRead(u,&U);CHKERRQ(ierr); 3891 ierr = PetscDrawAxisGetLimits(ictx->axis,&xl,&xr,&yl,&yr);CHKERRQ(ierr); 3892 if ((PetscRealPart(U[0]) < xl) || (PetscRealPart(U[1]) < yl) || (PetscRealPart(U[0]) > xr) || (PetscRealPart(U[1]) > yr)) { 3893 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3894 PetscFunctionReturn(0); 3895 } 3896 if (!step) ictx->color++; 3897 ierr = PetscDrawPoint(draw,PetscRealPart(U[0]),PetscRealPart(U[1]),ictx->color);CHKERRQ(ierr); 3898 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3899 3900 if (ictx->showtimestepandtime) { 3901 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3902 ierr = PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);CHKERRQ(ierr); 3903 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3904 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3905 w = xl + .5*(xr - xl) - .5*len*tw; 3906 h = yl + .95*(yr - yl); 3907 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3908 } 3909 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3910 PetscFunctionReturn(0); 3911 } 3912 3913 3914 #undef __FUNCT__ 3915 #define __FUNCT__ "TSMonitorDrawCtxDestroy" 3916 /*@C 3917 TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution() 3918 3919 Collective on TS 3920 3921 Input Parameters: 3922 . ctx - the monitor context 3923 3924 Level: intermediate 3925 3926 .keywords: TS, vector, monitor, view 3927 3928 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError() 3929 @*/ 3930 PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx) 3931 { 3932 PetscErrorCode ierr; 3933 3934 PetscFunctionBegin; 3935 ierr = PetscDrawAxisDestroy(&(*ictx)->axis);CHKERRQ(ierr); 3936 ierr = PetscViewerDestroy(&(*ictx)->viewer);CHKERRQ(ierr); 3937 ierr = VecDestroy(&(*ictx)->initialsolution);CHKERRQ(ierr); 3938 ierr = PetscFree(*ictx);CHKERRQ(ierr); 3939 PetscFunctionReturn(0); 3940 } 3941 3942 #undef __FUNCT__ 3943 #define __FUNCT__ "TSMonitorDrawCtxCreate" 3944 /*@C 3945 TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx 3946 3947 Collective on TS 3948 3949 Input Parameter: 3950 . ts - time-step context 3951 3952 Output Patameter: 3953 . ctx - the monitor context 3954 3955 Options Database: 3956 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3957 3958 Level: intermediate 3959 3960 .keywords: TS, vector, monitor, view 3961 3962 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx() 3963 @*/ 3964 PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx) 3965 { 3966 PetscErrorCode ierr; 3967 3968 PetscFunctionBegin; 3969 ierr = PetscNew(ctx);CHKERRQ(ierr); 3970 ierr = PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);CHKERRQ(ierr); 3971 ierr = PetscViewerSetFromOptions((*ctx)->viewer);CHKERRQ(ierr); 3972 3973 (*ctx)->howoften = howoften; 3974 (*ctx)->showinitial = PETSC_FALSE; 3975 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);CHKERRQ(ierr); 3976 3977 (*ctx)->showtimestepandtime = PETSC_FALSE; 3978 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);CHKERRQ(ierr); 3979 (*ctx)->color = PETSC_DRAW_WHITE; 3980 PetscFunctionReturn(0); 3981 } 3982 3983 #undef __FUNCT__ 3984 #define __FUNCT__ "TSMonitorDrawError" 3985 /*@C 3986 TSMonitorDrawError - Monitors progress of the TS solvers by calling 3987 VecView() for the error at each timestep 3988 3989 Collective on TS 3990 3991 Input Parameters: 3992 + ts - the TS context 3993 . step - current time-step 3994 . ptime - current time 3995 - dummy - either a viewer or NULL 3996 3997 Level: intermediate 3998 3999 .keywords: TS, vector, monitor, view 4000 4001 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4002 @*/ 4003 PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4004 { 4005 PetscErrorCode ierr; 4006 TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy; 4007 PetscViewer viewer = ctx->viewer; 4008 Vec work; 4009 4010 PetscFunctionBegin; 4011 if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 4012 ierr = VecDuplicate(u,&work);CHKERRQ(ierr); 4013 ierr = TSComputeSolutionFunction(ts,ptime,work);CHKERRQ(ierr); 4014 ierr = VecAXPY(work,-1.0,u);CHKERRQ(ierr); 4015 ierr = VecView(work,viewer);CHKERRQ(ierr); 4016 ierr = VecDestroy(&work);CHKERRQ(ierr); 4017 PetscFunctionReturn(0); 4018 } 4019 4020 #include <petsc-private/dmimpl.h> 4021 #undef __FUNCT__ 4022 #define __FUNCT__ "TSSetDM" 4023 /*@ 4024 TSSetDM - Sets the DM that may be used by some preconditioners 4025 4026 Logically Collective on TS and DM 4027 4028 Input Parameters: 4029 + ts - the preconditioner context 4030 - dm - the dm 4031 4032 Level: intermediate 4033 4034 4035 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 4036 @*/ 4037 PetscErrorCode TSSetDM(TS ts,DM dm) 4038 { 4039 PetscErrorCode ierr; 4040 SNES snes; 4041 DMTS tsdm; 4042 4043 PetscFunctionBegin; 4044 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4045 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 4046 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 4047 if (ts->dm->dmts && !dm->dmts) { 4048 ierr = DMCopyDMTS(ts->dm,dm);CHKERRQ(ierr); 4049 ierr = DMGetDMTS(ts->dm,&tsdm);CHKERRQ(ierr); 4050 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 4051 tsdm->originaldm = dm; 4052 } 4053 } 4054 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 4055 } 4056 ts->dm = dm; 4057 4058 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4059 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 4060 PetscFunctionReturn(0); 4061 } 4062 4063 #undef __FUNCT__ 4064 #define __FUNCT__ "TSGetDM" 4065 /*@ 4066 TSGetDM - Gets the DM that may be used by some preconditioners 4067 4068 Not Collective 4069 4070 Input Parameter: 4071 . ts - the preconditioner context 4072 4073 Output Parameter: 4074 . dm - the dm 4075 4076 Level: intermediate 4077 4078 4079 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 4080 @*/ 4081 PetscErrorCode TSGetDM(TS ts,DM *dm) 4082 { 4083 PetscErrorCode ierr; 4084 4085 PetscFunctionBegin; 4086 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4087 if (!ts->dm) { 4088 ierr = DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);CHKERRQ(ierr); 4089 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 4090 } 4091 *dm = ts->dm; 4092 PetscFunctionReturn(0); 4093 } 4094 4095 #undef __FUNCT__ 4096 #define __FUNCT__ "SNESTSFormFunction" 4097 /*@ 4098 SNESTSFormFunction - Function to evaluate nonlinear residual 4099 4100 Logically Collective on SNES 4101 4102 Input Parameter: 4103 + snes - nonlinear solver 4104 . U - the current state at which to evaluate the residual 4105 - ctx - user context, must be a TS 4106 4107 Output Parameter: 4108 . F - the nonlinear residual 4109 4110 Notes: 4111 This function is not normally called by users and is automatically registered with the SNES used by TS. 4112 It is most frequently passed to MatFDColoringSetFunction(). 4113 4114 Level: advanced 4115 4116 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 4117 @*/ 4118 PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx) 4119 { 4120 TS ts = (TS)ctx; 4121 PetscErrorCode ierr; 4122 4123 PetscFunctionBegin; 4124 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4125 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 4126 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 4127 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 4128 ierr = (ts->ops->snesfunction)(snes,U,F,ts);CHKERRQ(ierr); 4129 PetscFunctionReturn(0); 4130 } 4131 4132 #undef __FUNCT__ 4133 #define __FUNCT__ "SNESTSFormJacobian" 4134 /*@ 4135 SNESTSFormJacobian - Function to evaluate the Jacobian 4136 4137 Collective on SNES 4138 4139 Input Parameter: 4140 + snes - nonlinear solver 4141 . U - the current state at which to evaluate the residual 4142 - ctx - user context, must be a TS 4143 4144 Output Parameter: 4145 + A - the Jacobian 4146 . B - the preconditioning matrix (may be the same as A) 4147 - flag - indicates any structure change in the matrix 4148 4149 Notes: 4150 This function is not normally called by users and is automatically registered with the SNES used by TS. 4151 4152 Level: developer 4153 4154 .seealso: SNESSetJacobian() 4155 @*/ 4156 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat A,Mat B,void *ctx) 4157 { 4158 TS ts = (TS)ctx; 4159 PetscErrorCode ierr; 4160 4161 PetscFunctionBegin; 4162 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4163 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 4164 PetscValidPointer(A,3); 4165 PetscValidHeaderSpecific(A,MAT_CLASSID,3); 4166 PetscValidPointer(B,4); 4167 PetscValidHeaderSpecific(B,MAT_CLASSID,4); 4168 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 4169 ierr = (ts->ops->snesjacobian)(snes,U,A,B,ts);CHKERRQ(ierr); 4170 PetscFunctionReturn(0); 4171 } 4172 4173 #undef __FUNCT__ 4174 #define __FUNCT__ "TSComputeRHSFunctionLinear" 4175 /*@C 4176 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 4177 4178 Collective on TS 4179 4180 Input Arguments: 4181 + ts - time stepping context 4182 . t - time at which to evaluate 4183 . U - state at which to evaluate 4184 - ctx - context 4185 4186 Output Arguments: 4187 . F - right hand side 4188 4189 Level: intermediate 4190 4191 Notes: 4192 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 4193 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 4194 4195 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 4196 @*/ 4197 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 4198 { 4199 PetscErrorCode ierr; 4200 Mat Arhs,Brhs; 4201 4202 PetscFunctionBegin; 4203 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 4204 ierr = TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);CHKERRQ(ierr); 4205 ierr = MatMult(Arhs,U,F);CHKERRQ(ierr); 4206 PetscFunctionReturn(0); 4207 } 4208 4209 #undef __FUNCT__ 4210 #define __FUNCT__ "TSComputeRHSJacobianConstant" 4211 /*@C 4212 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 4213 4214 Collective on TS 4215 4216 Input Arguments: 4217 + ts - time stepping context 4218 . t - time at which to evaluate 4219 . U - state at which to evaluate 4220 - ctx - context 4221 4222 Output Arguments: 4223 + A - pointer to operator 4224 . B - pointer to preconditioning matrix 4225 - flg - matrix structure flag 4226 4227 Level: intermediate 4228 4229 Notes: 4230 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 4231 4232 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 4233 @*/ 4234 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx) 4235 { 4236 PetscFunctionBegin; 4237 PetscFunctionReturn(0); 4238 } 4239 4240 #undef __FUNCT__ 4241 #define __FUNCT__ "TSComputeIFunctionLinear" 4242 /*@C 4243 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 4244 4245 Collective on TS 4246 4247 Input Arguments: 4248 + ts - time stepping context 4249 . t - time at which to evaluate 4250 . U - state at which to evaluate 4251 . Udot - time derivative of state vector 4252 - ctx - context 4253 4254 Output Arguments: 4255 . F - left hand side 4256 4257 Level: intermediate 4258 4259 Notes: 4260 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 4261 user is required to write their own TSComputeIFunction. 4262 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 4263 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 4264 4265 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 4266 @*/ 4267 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 4268 { 4269 PetscErrorCode ierr; 4270 Mat A,B; 4271 4272 PetscFunctionBegin; 4273 ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr); 4274 ierr = TSComputeIJacobian(ts,t,U,Udot,1.0,A,B,PETSC_TRUE);CHKERRQ(ierr); 4275 ierr = MatMult(A,Udot,F);CHKERRQ(ierr); 4276 PetscFunctionReturn(0); 4277 } 4278 4279 #undef __FUNCT__ 4280 #define __FUNCT__ "TSComputeIJacobianConstant" 4281 /*@C 4282 TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE 4283 4284 Collective on TS 4285 4286 Input Arguments: 4287 + ts - time stepping context 4288 . t - time at which to evaluate 4289 . U - state at which to evaluate 4290 . Udot - time derivative of state vector 4291 . shift - shift to apply 4292 - ctx - context 4293 4294 Output Arguments: 4295 + A - pointer to operator 4296 . B - pointer to preconditioning matrix 4297 - flg - matrix structure flag 4298 4299 Level: advanced 4300 4301 Notes: 4302 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 4303 4304 It is only appropriate for problems of the form 4305 4306 $ M Udot = F(U,t) 4307 4308 where M is constant and F is non-stiff. The user must pass M to TSSetIJacobian(). The current implementation only 4309 works with IMEX time integration methods such as TSROSW and TSARKIMEX, since there is no support for de-constructing 4310 an implicit operator of the form 4311 4312 $ shift*M + J 4313 4314 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 4315 a copy of M or reassemble it when requested. 4316 4317 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 4318 @*/ 4319 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,void *ctx) 4320 { 4321 PetscErrorCode ierr; 4322 4323 PetscFunctionBegin; 4324 ierr = MatScale(A, shift / ts->ijacobian.shift);CHKERRQ(ierr); 4325 ts->ijacobian.shift = shift; 4326 PetscFunctionReturn(0); 4327 } 4328 4329 #undef __FUNCT__ 4330 #define __FUNCT__ "TSGetEquationType" 4331 /*@ 4332 TSGetEquationType - Gets the type of the equation that TS is solving. 4333 4334 Not Collective 4335 4336 Input Parameter: 4337 . ts - the TS context 4338 4339 Output Parameter: 4340 . equation_type - see TSEquationType 4341 4342 Level: beginner 4343 4344 .keywords: TS, equation type 4345 4346 .seealso: TSSetEquationType(), TSEquationType 4347 @*/ 4348 PetscErrorCode TSGetEquationType(TS ts,TSEquationType *equation_type) 4349 { 4350 PetscFunctionBegin; 4351 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4352 PetscValidPointer(equation_type,2); 4353 *equation_type = ts->equation_type; 4354 PetscFunctionReturn(0); 4355 } 4356 4357 #undef __FUNCT__ 4358 #define __FUNCT__ "TSSetEquationType" 4359 /*@ 4360 TSSetEquationType - Sets the type of the equation that TS is solving. 4361 4362 Not Collective 4363 4364 Input Parameter: 4365 + ts - the TS context 4366 . equation_type - see TSEquationType 4367 4368 Level: advanced 4369 4370 .keywords: TS, equation type 4371 4372 .seealso: TSGetEquationType(), TSEquationType 4373 @*/ 4374 PetscErrorCode TSSetEquationType(TS ts,TSEquationType equation_type) 4375 { 4376 PetscFunctionBegin; 4377 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4378 ts->equation_type = equation_type; 4379 PetscFunctionReturn(0); 4380 } 4381 4382 #undef __FUNCT__ 4383 #define __FUNCT__ "TSGetConvergedReason" 4384 /*@ 4385 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 4386 4387 Not Collective 4388 4389 Input Parameter: 4390 . ts - the TS context 4391 4392 Output Parameter: 4393 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 4394 manual pages for the individual convergence tests for complete lists 4395 4396 Level: beginner 4397 4398 Notes: 4399 Can only be called after the call to TSSolve() is complete. 4400 4401 .keywords: TS, nonlinear, set, convergence, test 4402 4403 .seealso: TSSetConvergenceTest(), TSConvergedReason 4404 @*/ 4405 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 4406 { 4407 PetscFunctionBegin; 4408 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4409 PetscValidPointer(reason,2); 4410 *reason = ts->reason; 4411 PetscFunctionReturn(0); 4412 } 4413 4414 #undef __FUNCT__ 4415 #define __FUNCT__ "TSSetConvergedReason" 4416 /*@ 4417 TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve. 4418 4419 Not Collective 4420 4421 Input Parameter: 4422 + ts - the TS context 4423 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 4424 manual pages for the individual convergence tests for complete lists 4425 4426 Level: advanced 4427 4428 Notes: 4429 Can only be called during TSSolve() is active. 4430 4431 .keywords: TS, nonlinear, set, convergence, test 4432 4433 .seealso: TSConvergedReason 4434 @*/ 4435 PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason) 4436 { 4437 PetscFunctionBegin; 4438 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4439 ts->reason = reason; 4440 PetscFunctionReturn(0); 4441 } 4442 4443 #undef __FUNCT__ 4444 #define __FUNCT__ "TSGetSolveTime" 4445 /*@ 4446 TSGetSolveTime - Gets the time after a call to TSSolve() 4447 4448 Not Collective 4449 4450 Input Parameter: 4451 . ts - the TS context 4452 4453 Output Parameter: 4454 . ftime - the final time. This time should correspond to the final time set with TSSetDuration() 4455 4456 Level: beginner 4457 4458 Notes: 4459 Can only be called after the call to TSSolve() is complete. 4460 4461 .keywords: TS, nonlinear, set, convergence, test 4462 4463 .seealso: TSSetConvergenceTest(), TSConvergedReason 4464 @*/ 4465 PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime) 4466 { 4467 PetscFunctionBegin; 4468 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4469 PetscValidPointer(ftime,2); 4470 *ftime = ts->solvetime; 4471 PetscFunctionReturn(0); 4472 } 4473 4474 #undef __FUNCT__ 4475 #define __FUNCT__ "TSGetTotalSteps" 4476 /*@ 4477 TSGetTotalSteps - Gets the total number of steps done since the last call to TSSetUp() or TSCreate() 4478 4479 Not Collective 4480 4481 Input Parameter: 4482 . ts - the TS context 4483 4484 Output Parameter: 4485 . steps - the number of steps 4486 4487 Level: beginner 4488 4489 Notes: 4490 Includes the number of steps for all calls to TSSolve() since TSSetUp() was called 4491 4492 .keywords: TS, nonlinear, set, convergence, test 4493 4494 .seealso: TSSetConvergenceTest(), TSConvergedReason 4495 @*/ 4496 PetscErrorCode TSGetTotalSteps(TS ts,PetscInt *steps) 4497 { 4498 PetscFunctionBegin; 4499 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4500 PetscValidPointer(steps,2); 4501 *steps = ts->total_steps; 4502 PetscFunctionReturn(0); 4503 } 4504 4505 #undef __FUNCT__ 4506 #define __FUNCT__ "TSGetSNESIterations" 4507 /*@ 4508 TSGetSNESIterations - Gets the total number of nonlinear iterations 4509 used by the time integrator. 4510 4511 Not Collective 4512 4513 Input Parameter: 4514 . ts - TS context 4515 4516 Output Parameter: 4517 . nits - number of nonlinear iterations 4518 4519 Notes: 4520 This counter is reset to zero for each successive call to TSSolve(). 4521 4522 Level: intermediate 4523 4524 .keywords: TS, get, number, nonlinear, iterations 4525 4526 .seealso: TSGetKSPIterations() 4527 @*/ 4528 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 4529 { 4530 PetscFunctionBegin; 4531 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4532 PetscValidIntPointer(nits,2); 4533 *nits = ts->snes_its; 4534 PetscFunctionReturn(0); 4535 } 4536 4537 #undef __FUNCT__ 4538 #define __FUNCT__ "TSGetKSPIterations" 4539 /*@ 4540 TSGetKSPIterations - Gets the total number of linear iterations 4541 used by the time integrator. 4542 4543 Not Collective 4544 4545 Input Parameter: 4546 . ts - TS context 4547 4548 Output Parameter: 4549 . lits - number of linear iterations 4550 4551 Notes: 4552 This counter is reset to zero for each successive call to TSSolve(). 4553 4554 Level: intermediate 4555 4556 .keywords: TS, get, number, linear, iterations 4557 4558 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 4559 @*/ 4560 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 4561 { 4562 PetscFunctionBegin; 4563 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4564 PetscValidIntPointer(lits,2); 4565 *lits = ts->ksp_its; 4566 PetscFunctionReturn(0); 4567 } 4568 4569 #undef __FUNCT__ 4570 #define __FUNCT__ "TSGetStepRejections" 4571 /*@ 4572 TSGetStepRejections - Gets the total number of rejected steps. 4573 4574 Not Collective 4575 4576 Input Parameter: 4577 . ts - TS context 4578 4579 Output Parameter: 4580 . rejects - number of steps rejected 4581 4582 Notes: 4583 This counter is reset to zero for each successive call to TSSolve(). 4584 4585 Level: intermediate 4586 4587 .keywords: TS, get, number 4588 4589 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 4590 @*/ 4591 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 4592 { 4593 PetscFunctionBegin; 4594 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4595 PetscValidIntPointer(rejects,2); 4596 *rejects = ts->reject; 4597 PetscFunctionReturn(0); 4598 } 4599 4600 #undef __FUNCT__ 4601 #define __FUNCT__ "TSGetSNESFailures" 4602 /*@ 4603 TSGetSNESFailures - Gets the total number of failed SNES solves 4604 4605 Not Collective 4606 4607 Input Parameter: 4608 . ts - TS context 4609 4610 Output Parameter: 4611 . fails - number of failed nonlinear solves 4612 4613 Notes: 4614 This counter is reset to zero for each successive call to TSSolve(). 4615 4616 Level: intermediate 4617 4618 .keywords: TS, get, number 4619 4620 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 4621 @*/ 4622 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 4623 { 4624 PetscFunctionBegin; 4625 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4626 PetscValidIntPointer(fails,2); 4627 *fails = ts->num_snes_failures; 4628 PetscFunctionReturn(0); 4629 } 4630 4631 #undef __FUNCT__ 4632 #define __FUNCT__ "TSSetMaxStepRejections" 4633 /*@ 4634 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 4635 4636 Not Collective 4637 4638 Input Parameter: 4639 + ts - TS context 4640 - rejects - maximum number of rejected steps, pass -1 for unlimited 4641 4642 Notes: 4643 The counter is reset to zero for each step 4644 4645 Options Database Key: 4646 . -ts_max_reject - Maximum number of step rejections before a step fails 4647 4648 Level: intermediate 4649 4650 .keywords: TS, set, maximum, number 4651 4652 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 4653 @*/ 4654 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 4655 { 4656 PetscFunctionBegin; 4657 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4658 ts->max_reject = rejects; 4659 PetscFunctionReturn(0); 4660 } 4661 4662 #undef __FUNCT__ 4663 #define __FUNCT__ "TSSetMaxSNESFailures" 4664 /*@ 4665 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 4666 4667 Not Collective 4668 4669 Input Parameter: 4670 + ts - TS context 4671 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 4672 4673 Notes: 4674 The counter is reset to zero for each successive call to TSSolve(). 4675 4676 Options Database Key: 4677 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 4678 4679 Level: intermediate 4680 4681 .keywords: TS, set, maximum, number 4682 4683 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 4684 @*/ 4685 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 4686 { 4687 PetscFunctionBegin; 4688 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4689 ts->max_snes_failures = fails; 4690 PetscFunctionReturn(0); 4691 } 4692 4693 #undef __FUNCT__ 4694 #define __FUNCT__ "TSSetErrorIfStepFails" 4695 /*@ 4696 TSSetErrorIfStepFails - Error if no step succeeds 4697 4698 Not Collective 4699 4700 Input Parameter: 4701 + ts - TS context 4702 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 4703 4704 Options Database Key: 4705 . -ts_error_if_step_fails - Error if no step succeeds 4706 4707 Level: intermediate 4708 4709 .keywords: TS, set, error 4710 4711 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 4712 @*/ 4713 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 4714 { 4715 PetscFunctionBegin; 4716 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4717 ts->errorifstepfailed = err; 4718 PetscFunctionReturn(0); 4719 } 4720 4721 #undef __FUNCT__ 4722 #define __FUNCT__ "TSMonitorSolutionBinary" 4723 /*@C 4724 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 4725 4726 Collective on TS 4727 4728 Input Parameters: 4729 + ts - the TS context 4730 . step - current time-step 4731 . ptime - current time 4732 . u - current state 4733 - viewer - binary viewer 4734 4735 Level: intermediate 4736 4737 .keywords: TS, vector, monitor, view 4738 4739 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4740 @*/ 4741 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec u,void *viewer) 4742 { 4743 PetscErrorCode ierr; 4744 PetscViewer v = (PetscViewer)viewer; 4745 4746 PetscFunctionBegin; 4747 ierr = VecView(u,v);CHKERRQ(ierr); 4748 PetscFunctionReturn(0); 4749 } 4750 4751 #undef __FUNCT__ 4752 #define __FUNCT__ "TSMonitorSolutionVTK" 4753 /*@C 4754 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 4755 4756 Collective on TS 4757 4758 Input Parameters: 4759 + ts - the TS context 4760 . step - current time-step 4761 . ptime - current time 4762 . u - current state 4763 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4764 4765 Level: intermediate 4766 4767 Notes: 4768 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. 4769 These are named according to the file name template. 4770 4771 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 4772 4773 .keywords: TS, vector, monitor, view 4774 4775 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4776 @*/ 4777 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate) 4778 { 4779 PetscErrorCode ierr; 4780 char filename[PETSC_MAX_PATH_LEN]; 4781 PetscViewer viewer; 4782 4783 PetscFunctionBegin; 4784 ierr = PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);CHKERRQ(ierr); 4785 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 4786 ierr = VecView(u,viewer);CHKERRQ(ierr); 4787 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4788 PetscFunctionReturn(0); 4789 } 4790 4791 #undef __FUNCT__ 4792 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 4793 /*@C 4794 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 4795 4796 Collective on TS 4797 4798 Input Parameters: 4799 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4800 4801 Level: intermediate 4802 4803 Note: 4804 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 4805 4806 .keywords: TS, vector, monitor, view 4807 4808 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 4809 @*/ 4810 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 4811 { 4812 PetscErrorCode ierr; 4813 4814 PetscFunctionBegin; 4815 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 4816 PetscFunctionReturn(0); 4817 } 4818 4819 #undef __FUNCT__ 4820 #define __FUNCT__ "TSGetAdapt" 4821 /*@ 4822 TSGetAdapt - Get the adaptive controller context for the current method 4823 4824 Collective on TS if controller has not been created yet 4825 4826 Input Arguments: 4827 . ts - time stepping context 4828 4829 Output Arguments: 4830 . adapt - adaptive controller 4831 4832 Level: intermediate 4833 4834 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 4835 @*/ 4836 PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt) 4837 { 4838 PetscErrorCode ierr; 4839 4840 PetscFunctionBegin; 4841 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4842 PetscValidPointer(adapt,2); 4843 if (!ts->adapt) { 4844 ierr = TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);CHKERRQ(ierr); 4845 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->adapt);CHKERRQ(ierr); 4846 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 4847 } 4848 *adapt = ts->adapt; 4849 PetscFunctionReturn(0); 4850 } 4851 4852 #undef __FUNCT__ 4853 #define __FUNCT__ "TSSetTolerances" 4854 /*@ 4855 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 4856 4857 Logically Collective 4858 4859 Input Arguments: 4860 + ts - time integration context 4861 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 4862 . vatol - vector of absolute tolerances or NULL, used in preference to atol if present 4863 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 4864 - vrtol - vector of relative tolerances or NULL, used in preference to atol if present 4865 4866 Options Database keys: 4867 + -ts_rtol <rtol> - relative tolerance for local truncation error 4868 - -ts_atol <atol> Absolute tolerance for local truncation error 4869 4870 Level: beginner 4871 4872 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 4873 @*/ 4874 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 4875 { 4876 PetscErrorCode ierr; 4877 4878 PetscFunctionBegin; 4879 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 4880 if (vatol) { 4881 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 4882 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 4883 4884 ts->vatol = vatol; 4885 } 4886 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 4887 if (vrtol) { 4888 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 4889 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 4890 4891 ts->vrtol = vrtol; 4892 } 4893 PetscFunctionReturn(0); 4894 } 4895 4896 #undef __FUNCT__ 4897 #define __FUNCT__ "TSGetTolerances" 4898 /*@ 4899 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 4900 4901 Logically Collective 4902 4903 Input Arguments: 4904 . ts - time integration context 4905 4906 Output Arguments: 4907 + atol - scalar absolute tolerances, NULL to ignore 4908 . vatol - vector of absolute tolerances, NULL to ignore 4909 . rtol - scalar relative tolerances, NULL to ignore 4910 - vrtol - vector of relative tolerances, NULL to ignore 4911 4912 Level: beginner 4913 4914 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 4915 @*/ 4916 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 4917 { 4918 PetscFunctionBegin; 4919 if (atol) *atol = ts->atol; 4920 if (vatol) *vatol = ts->vatol; 4921 if (rtol) *rtol = ts->rtol; 4922 if (vrtol) *vrtol = ts->vrtol; 4923 PetscFunctionReturn(0); 4924 } 4925 4926 #undef __FUNCT__ 4927 #define __FUNCT__ "TSErrorNormWRMS" 4928 /*@ 4929 TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state 4930 4931 Collective on TS 4932 4933 Input Arguments: 4934 + ts - time stepping context 4935 - Y - state vector to be compared to ts->vec_sol 4936 4937 Output Arguments: 4938 . norm - weighted norm, a value of 1.0 is considered small 4939 4940 Level: developer 4941 4942 .seealso: TSSetTolerances() 4943 @*/ 4944 PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm) 4945 { 4946 PetscErrorCode ierr; 4947 PetscInt i,n,N; 4948 const PetscScalar *u,*y; 4949 Vec U; 4950 PetscReal sum,gsum; 4951 4952 PetscFunctionBegin; 4953 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4954 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 4955 PetscValidPointer(norm,3); 4956 U = ts->vec_sol; 4957 PetscCheckSameTypeAndComm(U,1,Y,2); 4958 if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 4959 4960 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 4961 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 4962 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 4963 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 4964 sum = 0.; 4965 if (ts->vatol && ts->vrtol) { 4966 const PetscScalar *atol,*rtol; 4967 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4968 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4969 for (i=0; i<n; i++) { 4970 PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4971 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4972 } 4973 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4974 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4975 } else if (ts->vatol) { /* vector atol, scalar rtol */ 4976 const PetscScalar *atol; 4977 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4978 for (i=0; i<n; i++) { 4979 PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4980 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4981 } 4982 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4983 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 4984 const PetscScalar *rtol; 4985 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4986 for (i=0; i<n; i++) { 4987 PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4988 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4989 } 4990 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4991 } else { /* scalar atol, scalar rtol */ 4992 for (i=0; i<n; i++) { 4993 PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4994 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4995 } 4996 } 4997 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 4998 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 4999 5000 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 5001 *norm = PetscSqrtReal(gsum / N); 5002 if (PetscIsInfOrNanReal(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 5003 PetscFunctionReturn(0); 5004 } 5005 5006 #undef __FUNCT__ 5007 #define __FUNCT__ "TSSetCFLTimeLocal" 5008 /*@ 5009 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 5010 5011 Logically Collective on TS 5012 5013 Input Arguments: 5014 + ts - time stepping context 5015 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 5016 5017 Note: 5018 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 5019 5020 Level: intermediate 5021 5022 .seealso: TSGetCFLTime(), TSADAPTCFL 5023 @*/ 5024 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 5025 { 5026 PetscFunctionBegin; 5027 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5028 ts->cfltime_local = cfltime; 5029 ts->cfltime = -1.; 5030 PetscFunctionReturn(0); 5031 } 5032 5033 #undef __FUNCT__ 5034 #define __FUNCT__ "TSGetCFLTime" 5035 /*@ 5036 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 5037 5038 Collective on TS 5039 5040 Input Arguments: 5041 . ts - time stepping context 5042 5043 Output Arguments: 5044 . cfltime - maximum stable time step for forward Euler 5045 5046 Level: advanced 5047 5048 .seealso: TSSetCFLTimeLocal() 5049 @*/ 5050 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 5051 { 5052 PetscErrorCode ierr; 5053 5054 PetscFunctionBegin; 5055 if (ts->cfltime < 0) { 5056 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 5057 } 5058 *cfltime = ts->cfltime; 5059 PetscFunctionReturn(0); 5060 } 5061 5062 #undef __FUNCT__ 5063 #define __FUNCT__ "TSVISetVariableBounds" 5064 /*@ 5065 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 5066 5067 Input Parameters: 5068 . ts - the TS context. 5069 . xl - lower bound. 5070 . xu - upper bound. 5071 5072 Notes: 5073 If this routine is not called then the lower and upper bounds are set to 5074 PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp(). 5075 5076 Level: advanced 5077 5078 @*/ 5079 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 5080 { 5081 PetscErrorCode ierr; 5082 SNES snes; 5083 5084 PetscFunctionBegin; 5085 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 5086 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 5087 PetscFunctionReturn(0); 5088 } 5089 5090 #if defined(PETSC_HAVE_MATLAB_ENGINE) 5091 #include <mex.h> 5092 5093 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 5094 5095 #undef __FUNCT__ 5096 #define __FUNCT__ "TSComputeFunction_Matlab" 5097 /* 5098 TSComputeFunction_Matlab - Calls the function that has been set with 5099 TSSetFunctionMatlab(). 5100 5101 Collective on TS 5102 5103 Input Parameters: 5104 + snes - the TS context 5105 - u - input vector 5106 5107 Output Parameter: 5108 . y - function vector, as set by TSSetFunction() 5109 5110 Notes: 5111 TSComputeFunction() is typically used within nonlinear solvers 5112 implementations, so most users would not generally call this routine 5113 themselves. 5114 5115 Level: developer 5116 5117 .keywords: TS, nonlinear, compute, function 5118 5119 .seealso: TSSetFunction(), TSGetFunction() 5120 */ 5121 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx) 5122 { 5123 PetscErrorCode ierr; 5124 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5125 int nlhs = 1,nrhs = 7; 5126 mxArray *plhs[1],*prhs[7]; 5127 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 5128 5129 PetscFunctionBegin; 5130 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 5131 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 5132 PetscValidHeaderSpecific(udot,VEC_CLASSID,4); 5133 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 5134 PetscCheckSameComm(snes,1,u,3); 5135 PetscCheckSameComm(snes,1,y,5); 5136 5137 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5138 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5139 ierr = PetscMemcpy(&lxdot,&udot,sizeof(udot));CHKERRQ(ierr); 5140 ierr = PetscMemcpy(&ly,&y,sizeof(u));CHKERRQ(ierr); 5141 5142 prhs[0] = mxCreateDoubleScalar((double)ls); 5143 prhs[1] = mxCreateDoubleScalar(time); 5144 prhs[2] = mxCreateDoubleScalar((double)lx); 5145 prhs[3] = mxCreateDoubleScalar((double)lxdot); 5146 prhs[4] = mxCreateDoubleScalar((double)ly); 5147 prhs[5] = mxCreateString(sctx->funcname); 5148 prhs[6] = sctx->ctx; 5149 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 5150 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5151 mxDestroyArray(prhs[0]); 5152 mxDestroyArray(prhs[1]); 5153 mxDestroyArray(prhs[2]); 5154 mxDestroyArray(prhs[3]); 5155 mxDestroyArray(prhs[4]); 5156 mxDestroyArray(prhs[5]); 5157 mxDestroyArray(plhs[0]); 5158 PetscFunctionReturn(0); 5159 } 5160 5161 5162 #undef __FUNCT__ 5163 #define __FUNCT__ "TSSetFunctionMatlab" 5164 /* 5165 TSSetFunctionMatlab - Sets the function evaluation routine and function 5166 vector for use by the TS routines in solving ODEs 5167 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5168 5169 Logically Collective on TS 5170 5171 Input Parameters: 5172 + ts - the TS context 5173 - func - function evaluation routine 5174 5175 Calling sequence of func: 5176 $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx); 5177 5178 Level: beginner 5179 5180 .keywords: TS, nonlinear, set, function 5181 5182 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5183 */ 5184 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 5185 { 5186 PetscErrorCode ierr; 5187 TSMatlabContext *sctx; 5188 5189 PetscFunctionBegin; 5190 /* currently sctx is memory bleed */ 5191 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5192 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5193 /* 5194 This should work, but it doesn't 5195 sctx->ctx = ctx; 5196 mexMakeArrayPersistent(sctx->ctx); 5197 */ 5198 sctx->ctx = mxDuplicateArray(ctx); 5199 5200 ierr = TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 5201 PetscFunctionReturn(0); 5202 } 5203 5204 #undef __FUNCT__ 5205 #define __FUNCT__ "TSComputeJacobian_Matlab" 5206 /* 5207 TSComputeJacobian_Matlab - Calls the function that has been set with 5208 TSSetJacobianMatlab(). 5209 5210 Collective on TS 5211 5212 Input Parameters: 5213 + ts - the TS context 5214 . u - input vector 5215 . A, B - the matrices 5216 - ctx - user context 5217 5218 Level: developer 5219 5220 .keywords: TS, nonlinear, compute, function 5221 5222 .seealso: TSSetFunction(), TSGetFunction() 5223 @*/ 5224 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat A,Mat B,void *ctx) 5225 { 5226 PetscErrorCode ierr; 5227 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5228 int nlhs = 2,nrhs = 9; 5229 mxArray *plhs[2],*prhs[9]; 5230 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 5231 5232 PetscFunctionBegin; 5233 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5234 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 5235 5236 /* call Matlab function in ctx with arguments u and y */ 5237 5238 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 5239 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5240 ierr = PetscMemcpy(&lxdot,&udot,sizeof(u));CHKERRQ(ierr); 5241 ierr = PetscMemcpy(&lA,A,sizeof(u));CHKERRQ(ierr); 5242 ierr = PetscMemcpy(&lB,B,sizeof(u));CHKERRQ(ierr); 5243 5244 prhs[0] = mxCreateDoubleScalar((double)ls); 5245 prhs[1] = mxCreateDoubleScalar((double)time); 5246 prhs[2] = mxCreateDoubleScalar((double)lx); 5247 prhs[3] = mxCreateDoubleScalar((double)lxdot); 5248 prhs[4] = mxCreateDoubleScalar((double)shift); 5249 prhs[5] = mxCreateDoubleScalar((double)lA); 5250 prhs[6] = mxCreateDoubleScalar((double)lB); 5251 prhs[7] = mxCreateString(sctx->funcname); 5252 prhs[8] = sctx->ctx; 5253 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 5254 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5255 mxDestroyArray(prhs[0]); 5256 mxDestroyArray(prhs[1]); 5257 mxDestroyArray(prhs[2]); 5258 mxDestroyArray(prhs[3]); 5259 mxDestroyArray(prhs[4]); 5260 mxDestroyArray(prhs[5]); 5261 mxDestroyArray(prhs[6]); 5262 mxDestroyArray(prhs[7]); 5263 mxDestroyArray(plhs[0]); 5264 mxDestroyArray(plhs[1]); 5265 PetscFunctionReturn(0); 5266 } 5267 5268 5269 #undef __FUNCT__ 5270 #define __FUNCT__ "TSSetJacobianMatlab" 5271 /* 5272 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5273 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 5274 5275 Logically Collective on TS 5276 5277 Input Parameters: 5278 + ts - the TS context 5279 . A,B - Jacobian matrices 5280 . func - function evaluation routine 5281 - ctx - user context 5282 5283 Calling sequence of func: 5284 $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx); 5285 5286 5287 Level: developer 5288 5289 .keywords: TS, nonlinear, set, function 5290 5291 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5292 */ 5293 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 5294 { 5295 PetscErrorCode ierr; 5296 TSMatlabContext *sctx; 5297 5298 PetscFunctionBegin; 5299 /* currently sctx is memory bleed */ 5300 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5301 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5302 /* 5303 This should work, but it doesn't 5304 sctx->ctx = ctx; 5305 mexMakeArrayPersistent(sctx->ctx); 5306 */ 5307 sctx->ctx = mxDuplicateArray(ctx); 5308 5309 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5310 PetscFunctionReturn(0); 5311 } 5312 5313 #undef __FUNCT__ 5314 #define __FUNCT__ "TSMonitor_Matlab" 5315 /* 5316 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 5317 5318 Collective on TS 5319 5320 .seealso: TSSetFunction(), TSGetFunction() 5321 @*/ 5322 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx) 5323 { 5324 PetscErrorCode ierr; 5325 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5326 int nlhs = 1,nrhs = 6; 5327 mxArray *plhs[1],*prhs[6]; 5328 long long int lx = 0,ls = 0; 5329 5330 PetscFunctionBegin; 5331 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5332 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 5333 5334 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 5335 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5336 5337 prhs[0] = mxCreateDoubleScalar((double)ls); 5338 prhs[1] = mxCreateDoubleScalar((double)it); 5339 prhs[2] = mxCreateDoubleScalar((double)time); 5340 prhs[3] = mxCreateDoubleScalar((double)lx); 5341 prhs[4] = mxCreateString(sctx->funcname); 5342 prhs[5] = sctx->ctx; 5343 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 5344 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5345 mxDestroyArray(prhs[0]); 5346 mxDestroyArray(prhs[1]); 5347 mxDestroyArray(prhs[2]); 5348 mxDestroyArray(prhs[3]); 5349 mxDestroyArray(prhs[4]); 5350 mxDestroyArray(plhs[0]); 5351 PetscFunctionReturn(0); 5352 } 5353 5354 5355 #undef __FUNCT__ 5356 #define __FUNCT__ "TSMonitorSetMatlab" 5357 /* 5358 TSMonitorSetMatlab - Sets the monitor function from Matlab 5359 5360 Level: developer 5361 5362 .keywords: TS, nonlinear, set, function 5363 5364 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5365 */ 5366 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 5367 { 5368 PetscErrorCode ierr; 5369 TSMatlabContext *sctx; 5370 5371 PetscFunctionBegin; 5372 /* currently sctx is memory bleed */ 5373 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5374 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5375 /* 5376 This should work, but it doesn't 5377 sctx->ctx = ctx; 5378 mexMakeArrayPersistent(sctx->ctx); 5379 */ 5380 sctx->ctx = mxDuplicateArray(ctx); 5381 5382 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5383 PetscFunctionReturn(0); 5384 } 5385 #endif 5386 5387 #undef __FUNCT__ 5388 #define __FUNCT__ "TSMonitorLGSolution" 5389 /*@C 5390 TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector 5391 in a time based line graph 5392 5393 Collective on TS 5394 5395 Input Parameters: 5396 + ts - the TS context 5397 . step - current time-step 5398 . ptime - current time 5399 - lg - a line graph object 5400 5401 Options Database: 5402 . -ts_monitor_lg_solution_variables 5403 5404 Level: intermediate 5405 5406 Notes: each process in a parallel run displays its component solutions in a separate window 5407 5408 .keywords: TS, vector, monitor, view 5409 5410 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 5411 @*/ 5412 PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 5413 { 5414 PetscErrorCode ierr; 5415 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 5416 const PetscScalar *yy; 5417 PetscInt dim; 5418 Vec v; 5419 5420 PetscFunctionBegin; 5421 if (!step) { 5422 PetscDrawAxis axis; 5423 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 5424 ierr = PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");CHKERRQ(ierr); 5425 if (ctx->names && !ctx->displaynames) { 5426 char **displaynames; 5427 PetscBool flg; 5428 5429 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 5430 ierr = PetscMalloc((dim+1)*sizeof(char*),&displaynames);CHKERRQ(ierr); 5431 ierr = PetscMemzero(displaynames,(dim+1)*sizeof(char*));CHKERRQ(ierr); 5432 ierr = PetscOptionsGetStringArray(((PetscObject)ts)->prefix,"-ts_monitor_lg_solution_variables",displaynames,&dim,&flg);CHKERRQ(ierr); 5433 if (flg) { 5434 ierr = TSMonitorLGCtxSetDisplayVariables(ctx,(const char *const *)displaynames);CHKERRQ(ierr); 5435 } 5436 ierr = PetscStrArrayDestroy(&displaynames);CHKERRQ(ierr); 5437 } 5438 if (ctx->displaynames) { 5439 ierr = PetscDrawLGSetDimension(ctx->lg,ctx->ndisplayvariables);CHKERRQ(ierr); 5440 ierr = PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->displaynames);CHKERRQ(ierr); 5441 } else if (ctx->names) { 5442 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 5443 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 5444 ierr = PetscDrawLGSetLegend(ctx->lg,(const char *const *)ctx->names);CHKERRQ(ierr); 5445 } 5446 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 5447 } 5448 if (ctx->transform) { 5449 ierr = (*ctx->transform)(ctx->transformctx,u,&v);CHKERRQ(ierr); 5450 } else { 5451 v = u; 5452 } 5453 ierr = VecGetArrayRead(v,&yy);CHKERRQ(ierr); 5454 #if defined(PETSC_USE_COMPLEX) 5455 { 5456 PetscReal *yreal; 5457 PetscInt i,n; 5458 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 5459 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 5460 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 5461 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 5462 ierr = PetscFree(yreal);CHKERRQ(ierr); 5463 } 5464 #else 5465 if (ctx->displaynames) { 5466 PetscInt i; 5467 for (i=0; i<ctx->ndisplayvariables; i++) { 5468 ctx->displayvalues[i] = yy[ctx->displayvariables[i]]; 5469 } 5470 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,ctx->displayvalues);CHKERRQ(ierr); 5471 } else { 5472 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 5473 } 5474 #endif 5475 ierr = VecRestoreArrayRead(v,&yy);CHKERRQ(ierr); 5476 if (ctx->transform) { 5477 ierr = VecDestroy(&v);CHKERRQ(ierr); 5478 } 5479 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 5480 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 5481 } 5482 PetscFunctionReturn(0); 5483 } 5484 5485 5486 #undef __FUNCT__ 5487 #define __FUNCT__ "TSMonitorLGSetVariableNames" 5488 /*@C 5489 TSMonitorLGSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot 5490 5491 Collective on TS 5492 5493 Input Parameters: 5494 + ts - the TS context 5495 - names - the names of the components, final string must be NULL 5496 5497 Level: intermediate 5498 5499 .keywords: TS, vector, monitor, view 5500 5501 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGCtxSetVariableNames() 5502 @*/ 5503 PetscErrorCode TSMonitorLGSetVariableNames(TS ts,const char * const *names) 5504 { 5505 PetscErrorCode ierr; 5506 PetscInt i; 5507 5508 PetscFunctionBegin; 5509 for (i=0; i<ts->numbermonitors; i++) { 5510 if (ts->monitor[i] == TSMonitorLGSolution) { 5511 ierr = TSMonitorLGCtxSetVariableNames(ts->monitorcontext[i],names);CHKERRQ(ierr); 5512 break; 5513 } 5514 } 5515 PetscFunctionReturn(0); 5516 } 5517 5518 #undef __FUNCT__ 5519 #define __FUNCT__ "TSMonitorLGCtxSetVariableNames" 5520 /*@C 5521 TSMonitorLGCtxSetVariableNames - Sets the name of each component in the solution vector so that it may be displayed in the plot 5522 5523 Collective on TS 5524 5525 Input Parameters: 5526 + ts - the TS context 5527 - names - the names of the components, final string must be NULL 5528 5529 Level: intermediate 5530 5531 .keywords: TS, vector, monitor, view 5532 5533 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables(), TSMonitorLGSetVariableNames() 5534 @*/ 5535 PetscErrorCode TSMonitorLGCtxSetVariableNames(TSMonitorLGCtx ctx,const char * const *names) 5536 { 5537 PetscErrorCode ierr; 5538 5539 PetscFunctionBegin; 5540 ierr = PetscStrArrayDestroy(&ctx->names);CHKERRQ(ierr); 5541 ierr = PetscStrArrayallocpy(names,&ctx->names);CHKERRQ(ierr); 5542 PetscFunctionReturn(0); 5543 } 5544 5545 #undef __FUNCT__ 5546 #define __FUNCT__ "TSMonitorLGGetVariableNames" 5547 /*@C 5548 TSMonitorLGGetVariableNames - Gets the name of each component in the solution vector so that it may be displayed in the plot 5549 5550 Collective on TS 5551 5552 Input Parameter: 5553 . ts - the TS context 5554 5555 Output Parameter: 5556 . names - the names of the components, final string must be NULL 5557 5558 Level: intermediate 5559 5560 .keywords: TS, vector, monitor, view 5561 5562 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables() 5563 @*/ 5564 PetscErrorCode TSMonitorLGGetVariableNames(TS ts,const char *const **names) 5565 { 5566 PetscInt i; 5567 5568 PetscFunctionBegin; 5569 *names = NULL; 5570 for (i=0; i<ts->numbermonitors; i++) { 5571 if (ts->monitor[i] == TSMonitorLGSolution) { 5572 TSMonitorLGCtx ctx = ts->monitorcontext[i]; 5573 *names = (const char *const *)ctx->names; 5574 break; 5575 } 5576 } 5577 PetscFunctionReturn(0); 5578 } 5579 5580 #undef __FUNCT__ 5581 #define __FUNCT__ "TSMonitorLGCtxSetDisplayVariables" 5582 /*@C 5583 TSMonitorLGCtxSetDisplayVariables - Sets the variables that are to be display in the monitor 5584 5585 Collective on TS 5586 5587 Input Parameters: 5588 + ctx - the TSMonitorLG context 5589 . displaynames - the names of the components, final string must be NULL 5590 5591 Level: intermediate 5592 5593 .keywords: TS, vector, monitor, view 5594 5595 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames() 5596 @*/ 5597 PetscErrorCode TSMonitorLGCtxSetDisplayVariables(TSMonitorLGCtx ctx,const char * const *displaynames) 5598 { 5599 PetscInt j = 0,k; 5600 PetscErrorCode ierr; 5601 5602 PetscFunctionBegin; 5603 if (!ctx->names) PetscFunctionReturn(0); 5604 ierr = PetscStrArrayDestroy(&ctx->displaynames);CHKERRQ(ierr); 5605 ierr = PetscStrArrayallocpy(displaynames,&ctx->displaynames);CHKERRQ(ierr); 5606 while (displaynames[j]) j++; 5607 ctx->ndisplayvariables = j; 5608 ierr = PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvariables);CHKERRQ(ierr); 5609 ierr = PetscMalloc1(ctx->ndisplayvariables,&ctx->displayvalues);CHKERRQ(ierr); 5610 j = 0; 5611 while (displaynames[j]) { 5612 k = 0; 5613 while (ctx->names[k]) { 5614 PetscBool flg; 5615 ierr = PetscStrcmp(displaynames[j],ctx->names[k],&flg);CHKERRQ(ierr); 5616 if (flg) { 5617 ctx->displayvariables[j] = k; 5618 break; 5619 } 5620 k++; 5621 } 5622 j++; 5623 } 5624 PetscFunctionReturn(0); 5625 } 5626 5627 5628 #undef __FUNCT__ 5629 #define __FUNCT__ "TSMonitorLGSetDisplayVariables" 5630 /*@C 5631 TSMonitorLGSetDisplayVariables - Sets the variables that are to be display in the monitor 5632 5633 Collective on TS 5634 5635 Input Parameters: 5636 + ts - the TS context 5637 . displaynames - the names of the components, final string must be NULL 5638 5639 Level: intermediate 5640 5641 .keywords: TS, vector, monitor, view 5642 5643 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames() 5644 @*/ 5645 PetscErrorCode TSMonitorLGSetDisplayVariables(TS ts,const char * const *displaynames) 5646 { 5647 PetscInt i; 5648 PetscErrorCode ierr; 5649 5650 PetscFunctionBegin; 5651 for (i=0; i<ts->numbermonitors; i++) { 5652 if (ts->monitor[i] == TSMonitorLGSolution) { 5653 ierr = TSMonitorLGCtxSetDisplayVariables(ts->monitorcontext[i],displaynames);CHKERRQ(ierr); 5654 break; 5655 } 5656 } 5657 PetscFunctionReturn(0); 5658 } 5659 5660 #undef __FUNCT__ 5661 #define __FUNCT__ "TSMonitorLGSetTransform" 5662 /*@C 5663 TSMonitorLGSetTransform - Solution vector will be transformed by provided function before being displayed 5664 5665 Collective on TS 5666 5667 Input Parameters: 5668 + ts - the TS context 5669 . transform - the transform function 5670 . destroy - function to destroy the optional context 5671 - ctx - optional context used by transform function 5672 5673 Level: intermediate 5674 5675 .keywords: TS, vector, monitor, view 5676 5677 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGCtxSetTransform() 5678 @*/ 5679 PetscErrorCode TSMonitorLGSetTransform(TS ts,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx) 5680 { 5681 PetscInt i; 5682 PetscErrorCode ierr; 5683 5684 PetscFunctionBegin; 5685 for (i=0; i<ts->numbermonitors; i++) { 5686 if (ts->monitor[i] == TSMonitorLGSolution) { 5687 ierr = TSMonitorLGCtxSetTransform(ts->monitorcontext[i],transform,destroy,tctx);CHKERRQ(ierr); 5688 } 5689 } 5690 PetscFunctionReturn(0); 5691 } 5692 5693 #undef __FUNCT__ 5694 #define __FUNCT__ "TSMonitorLGCtxSetTransform" 5695 /*@C 5696 TSMonitorLGCtxSetTransform - Solution vector will be transformed by provided function before being displayed 5697 5698 Collective on TSLGCtx 5699 5700 Input Parameters: 5701 + ts - the TS context 5702 . transform - the transform function 5703 . destroy - function to destroy the optional context 5704 - ctx - optional context used by transform function 5705 5706 Level: intermediate 5707 5708 .keywords: TS, vector, monitor, view 5709 5710 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetVariableNames(), TSMonitorLGSetTransform() 5711 @*/ 5712 PetscErrorCode TSMonitorLGCtxSetTransform(TSMonitorLGCtx ctx,PetscErrorCode (*transform)(void*,Vec,Vec*),PetscErrorCode (*destroy)(void*),void *tctx) 5713 { 5714 PetscFunctionBegin; 5715 ctx->transform = transform; 5716 ctx->transformdestroy = destroy; 5717 ctx->transformctx = tctx; 5718 PetscFunctionReturn(0); 5719 } 5720 5721 #undef __FUNCT__ 5722 #define __FUNCT__ "TSMonitorLGError" 5723 /*@C 5724 TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector 5725 in a time based line graph 5726 5727 Collective on TS 5728 5729 Input Parameters: 5730 + ts - the TS context 5731 . step - current time-step 5732 . ptime - current time 5733 - lg - a line graph object 5734 5735 Level: intermediate 5736 5737 Notes: 5738 Only for sequential solves. 5739 5740 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 5741 5742 Options Database Keys: 5743 . -ts_monitor_lg_error - create a graphical monitor of error history 5744 5745 .keywords: TS, vector, monitor, view 5746 5747 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 5748 @*/ 5749 PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 5750 { 5751 PetscErrorCode ierr; 5752 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 5753 const PetscScalar *yy; 5754 Vec y; 5755 PetscInt dim; 5756 5757 PetscFunctionBegin; 5758 if (!step) { 5759 PetscDrawAxis axis; 5760 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 5761 ierr = PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");CHKERRQ(ierr); 5762 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 5763 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 5764 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 5765 } 5766 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 5767 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 5768 ierr = VecAXPY(y,-1.0,u);CHKERRQ(ierr); 5769 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 5770 #if defined(PETSC_USE_COMPLEX) 5771 { 5772 PetscReal *yreal; 5773 PetscInt i,n; 5774 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 5775 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 5776 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 5777 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 5778 ierr = PetscFree(yreal);CHKERRQ(ierr); 5779 } 5780 #else 5781 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 5782 #endif 5783 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 5784 ierr = VecDestroy(&y);CHKERRQ(ierr); 5785 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 5786 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 5787 } 5788 PetscFunctionReturn(0); 5789 } 5790 5791 #undef __FUNCT__ 5792 #define __FUNCT__ "TSMonitorLGSNESIterations" 5793 PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 5794 { 5795 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 5796 PetscReal x = ptime,y; 5797 PetscErrorCode ierr; 5798 PetscInt its; 5799 5800 PetscFunctionBegin; 5801 if (!n) { 5802 PetscDrawAxis axis; 5803 5804 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 5805 ierr = PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");CHKERRQ(ierr); 5806 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 5807 5808 ctx->snes_its = 0; 5809 } 5810 ierr = TSGetSNESIterations(ts,&its);CHKERRQ(ierr); 5811 y = its - ctx->snes_its; 5812 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 5813 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 5814 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 5815 } 5816 ctx->snes_its = its; 5817 PetscFunctionReturn(0); 5818 } 5819 5820 #undef __FUNCT__ 5821 #define __FUNCT__ "TSMonitorLGKSPIterations" 5822 PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 5823 { 5824 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 5825 PetscReal x = ptime,y; 5826 PetscErrorCode ierr; 5827 PetscInt its; 5828 5829 PetscFunctionBegin; 5830 if (!n) { 5831 PetscDrawAxis axis; 5832 5833 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 5834 ierr = PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");CHKERRQ(ierr); 5835 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 5836 5837 ctx->ksp_its = 0; 5838 } 5839 ierr = TSGetKSPIterations(ts,&its);CHKERRQ(ierr); 5840 y = its - ctx->ksp_its; 5841 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 5842 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 5843 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 5844 } 5845 ctx->ksp_its = its; 5846 PetscFunctionReturn(0); 5847 } 5848 5849 #undef __FUNCT__ 5850 #define __FUNCT__ "TSComputeLinearStability" 5851 /*@ 5852 TSComputeLinearStability - computes the linear stability function at a point 5853 5854 Collective on TS and Vec 5855 5856 Input Parameters: 5857 + ts - the TS context 5858 - xr,xi - real and imaginary part of input arguments 5859 5860 Output Parameters: 5861 . yr,yi - real and imaginary part of function value 5862 5863 Level: developer 5864 5865 .keywords: TS, compute 5866 5867 .seealso: TSSetRHSFunction(), TSComputeIFunction() 5868 @*/ 5869 PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 5870 { 5871 PetscErrorCode ierr; 5872 5873 PetscFunctionBegin; 5874 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5875 if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method"); 5876 ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr); 5877 PetscFunctionReturn(0); 5878 } 5879 5880 /* ------------------------------------------------------------------------*/ 5881 #undef __FUNCT__ 5882 #define __FUNCT__ "TSMonitorEnvelopeCtxCreate" 5883 /*@C 5884 TSMonitorEnvelopeCtxCreate - Creates a context for use with TSMonitorEnvelope() 5885 5886 Collective on TS 5887 5888 Input Parameters: 5889 . ts - the ODE solver object 5890 5891 Output Parameter: 5892 . ctx - the context 5893 5894 Level: intermediate 5895 5896 .keywords: TS, monitor, line graph, residual, seealso 5897 5898 .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 5899 5900 @*/ 5901 PetscErrorCode TSMonitorEnvelopeCtxCreate(TS ts,TSMonitorEnvelopeCtx *ctx) 5902 { 5903 PetscErrorCode ierr; 5904 5905 PetscFunctionBegin; 5906 ierr = PetscNew(ctx);CHKERRQ(ierr); 5907 PetscFunctionReturn(0); 5908 } 5909 5910 #undef __FUNCT__ 5911 #define __FUNCT__ "TSMonitorEnvelope" 5912 /*@C 5913 TSMonitorEnvelope - Monitors the maximum and minimum value of each component of the solution 5914 5915 Collective on TS 5916 5917 Input Parameters: 5918 + ts - the TS context 5919 . step - current time-step 5920 . ptime - current time 5921 - ctx - the envelope context 5922 5923 Options Database: 5924 . -ts_monitor_envelope 5925 5926 Level: intermediate 5927 5928 Notes: after a solve you can use TSMonitorEnvelopeGetBounds() to access the envelope 5929 5930 .keywords: TS, vector, monitor, view 5931 5932 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorEnvelopeGetBounds() 5933 @*/ 5934 PetscErrorCode TSMonitorEnvelope(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 5935 { 5936 PetscErrorCode ierr; 5937 TSMonitorEnvelopeCtx ctx = (TSMonitorEnvelopeCtx)dummy; 5938 5939 PetscFunctionBegin; 5940 if (!ctx->max) { 5941 ierr = VecDuplicate(u,&ctx->max);CHKERRQ(ierr); 5942 ierr = VecDuplicate(u,&ctx->min);CHKERRQ(ierr); 5943 ierr = VecCopy(u,ctx->max);CHKERRQ(ierr); 5944 ierr = VecCopy(u,ctx->min);CHKERRQ(ierr); 5945 } else { 5946 ierr = VecPointwiseMax(ctx->max,u,ctx->max);CHKERRQ(ierr); 5947 ierr = VecPointwiseMin(ctx->min,u,ctx->min);CHKERRQ(ierr); 5948 } 5949 PetscFunctionReturn(0); 5950 } 5951 5952 5953 #undef __FUNCT__ 5954 #define __FUNCT__ "TSMonitorEnvelopeGetBounds" 5955 /*@C 5956 TSMonitorEnvelopeGetBounds - Gets the bounds for the components of the solution 5957 5958 Collective on TS 5959 5960 Input Parameter: 5961 . ts - the TS context 5962 5963 Output Parameter: 5964 + max - the maximum values 5965 - min - the minimum values 5966 5967 Level: intermediate 5968 5969 .keywords: TS, vector, monitor, view 5970 5971 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorLGSetDisplayVariables() 5972 @*/ 5973 PetscErrorCode TSMonitorEnvelopeGetBounds(TS ts,Vec *max,Vec *min) 5974 { 5975 PetscInt i; 5976 5977 PetscFunctionBegin; 5978 if (max) *max = NULL; 5979 if (min) *min = NULL; 5980 for (i=0; i<ts->numbermonitors; i++) { 5981 if (ts->monitor[i] == TSMonitorEnvelope) { 5982 TSMonitorEnvelopeCtx ctx = ts->monitorcontext[i]; 5983 if (max) *max = ctx->max; 5984 if (min) *min = ctx->min; 5985 break; 5986 } 5987 } 5988 PetscFunctionReturn(0); 5989 } 5990 5991 #undef __FUNCT__ 5992 #define __FUNCT__ "TSMonitorEnvelopeCtxDestroy" 5993 /*@C 5994 TSMonitorEnvelopeCtxDestroy - Destroys a context that was created with TSMonitorEnvelopeCtxCreate(). 5995 5996 Collective on TSMonitorEnvelopeCtx 5997 5998 Input Parameter: 5999 . ctx - the monitor context 6000 6001 Level: intermediate 6002 6003 .keywords: TS, monitor, line graph, destroy 6004 6005 .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep(); 6006 @*/ 6007 PetscErrorCode TSMonitorEnvelopeCtxDestroy(TSMonitorEnvelopeCtx *ctx) 6008 { 6009 PetscErrorCode ierr; 6010 6011 PetscFunctionBegin; 6012 ierr = VecDestroy(&(*ctx)->min);CHKERRQ(ierr); 6013 ierr = VecDestroy(&(*ctx)->max);CHKERRQ(ierr); 6014 ierr = PetscFree(*ctx);CHKERRQ(ierr); 6015 PetscFunctionReturn(0); 6016 } 6017 6018 #undef __FUNCT__ 6019 #define __FUNCT__ "TSRollBack" 6020 /*@ 6021 TSRollBack - Rolls back one time step 6022 6023 Collective on TS 6024 6025 Input Parameter: 6026 . ts - the TS context obtained from TSCreate() 6027 6028 Level: advanced 6029 6030 .keywords: TS, timestep, rollback 6031 6032 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate() 6033 @*/ 6034 PetscErrorCode TSRollBack(TS ts) 6035 { 6036 PetscErrorCode ierr; 6037 6038 PetscFunctionBegin; 6039 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 6040 6041 if (!ts->ops->rollback) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRollBack not implemented for type '%s'",((PetscObject)ts)->type_name); 6042 ierr = (*ts->ops->rollback)(ts);CHKERRQ(ierr); 6043 ts->time_step = ts->ptime - ts->ptime_prev; 6044 ts->ptime = ts->ptime_prev; 6045 PetscFunctionReturn(0); 6046 } 6047 6048 #undef __FUNCT__ 6049 #define __FUNCT__ "TSGetStages" 6050 /*@ 6051 TSGetStages - Get the number of stages and stage values 6052 6053 Input Parameter: 6054 . ts - the TS context obtained from TSCreate() 6055 6056 Level: advanced 6057 6058 .keywords: TS, getstages 6059 6060 .seealso: TSCreate() 6061 @*/ 6062 PetscErrorCode TSGetStages(TS ts,PetscInt *ns, Vec **Y) 6063 { 6064 PetscErrorCode ierr; 6065 6066 PetscFunctionBegin; 6067 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 6068 PetscValidPointer(ns,2); 6069 6070 if (!ts->ops->getstages) *ns=0; 6071 else { 6072 ierr = (*ts->ops->getstages)(ts,ns,Y);CHKERRQ(ierr); 6073 } 6074 PetscFunctionReturn(0); 6075 } 6076 6077