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