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 PetscReal ptime_prev; 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 = TSSetUp(ts);CHKERRQ(ierr); 2625 2626 ts->reason = TS_CONVERGED_ITERATING; 2627 ts->ptime_prev = ts->ptime; 2628 2629 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2630 ierr = (*ts->ops->step)(ts);CHKERRQ(ierr); 2631 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2632 2633 ts->time_step_prev = ts->ptime - ts->ptime_prev; 2634 2635 if (ts->reason < 0) { 2636 if (ts->errorifstepfailed) { 2637 if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) { 2638 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]); 2639 } else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 2640 } 2641 } else if (!ts->reason) { 2642 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 2643 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 2644 } 2645 PetscFunctionReturn(0); 2646 } 2647 2648 #undef __FUNCT__ 2649 #define __FUNCT__ "TSEvaluateStep" 2650 /*@ 2651 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 2652 2653 Collective on TS 2654 2655 Input Arguments: 2656 + ts - time stepping context 2657 . order - desired order of accuracy 2658 - done - whether the step was evaluated at this order (pass NULL to generate an error if not available) 2659 2660 Output Arguments: 2661 . U - state at the end of the current step 2662 2663 Level: advanced 2664 2665 Notes: 2666 This function cannot be called until all stages have been evaluated. 2667 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. 2668 2669 .seealso: TSStep(), TSAdapt 2670 @*/ 2671 PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done) 2672 { 2673 PetscErrorCode ierr; 2674 2675 PetscFunctionBegin; 2676 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2677 PetscValidType(ts,1); 2678 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 2679 if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name); 2680 ierr = (*ts->ops->evaluatestep)(ts,order,U,done);CHKERRQ(ierr); 2681 PetscFunctionReturn(0); 2682 } 2683 2684 #undef __FUNCT__ 2685 #define __FUNCT__ "TSSolve" 2686 /*@ 2687 TSSolve - Steps the requested number of timesteps. 2688 2689 Collective on TS 2690 2691 Input Parameter: 2692 + ts - the TS context obtained from TSCreate() 2693 - u - the solution vector (can be null if TSSetSolution() was used, otherwise must contain the initial conditions) 2694 2695 Level: beginner 2696 2697 Notes: 2698 The final time returned by this function may be different from the time of the internally 2699 held state accessible by TSGetSolution() and TSGetTime() because the method may have 2700 stepped over the final time. 2701 2702 .keywords: TS, timestep, solve 2703 2704 .seealso: TSCreate(), TSSetSolution(), TSStep() 2705 @*/ 2706 PetscErrorCode TSSolve(TS ts,Vec u) 2707 { 2708 Vec solution; 2709 PetscErrorCode ierr; 2710 2711 PetscFunctionBegin; 2712 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2713 if (u) PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2714 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 */ 2715 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2716 if (!ts->vec_sol || u == ts->vec_sol) { 2717 ierr = VecDuplicate(u,&solution);CHKERRQ(ierr); 2718 ierr = TSSetSolution(ts,solution);CHKERRQ(ierr); 2719 ierr = VecDestroy(&solution);CHKERRQ(ierr); /* grant ownership */ 2720 } 2721 ierr = VecCopy(u,ts->vec_sol);CHKERRQ(ierr); 2722 } else if (u) { 2723 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 2724 } 2725 ierr = TSSetUp(ts);CHKERRQ(ierr); 2726 /* reset time step and iteration counters */ 2727 ts->steps = 0; 2728 ts->ksp_its = 0; 2729 ts->snes_its = 0; 2730 ts->num_snes_failures = 0; 2731 ts->reject = 0; 2732 ts->reason = TS_CONVERGED_ITERATING; 2733 2734 ierr = TSViewFromOptions(ts,NULL,"-ts_view_pre");CHKERRQ(ierr); 2735 2736 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 2737 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 2738 ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr); 2739 ts->solvetime = ts->ptime; 2740 } else { 2741 /* steps the requested number of timesteps. */ 2742 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 2743 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 2744 while (!ts->reason) { 2745 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2746 ierr = TSStep(ts);CHKERRQ(ierr); 2747 ierr = TSPostStep(ts);CHKERRQ(ierr); 2748 if (ts->event) { 2749 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 2750 } 2751 } 2752 if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) { 2753 ierr = TSInterpolate(ts,ts->max_time,u);CHKERRQ(ierr); 2754 ts->solvetime = ts->max_time; 2755 solution = u; 2756 } else { 2757 if (u) {ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr);} 2758 ts->solvetime = ts->ptime; 2759 solution = ts->vec_sol; 2760 } 2761 ierr = TSMonitor(ts,ts->steps,ts->solvetime,solution);CHKERRQ(ierr); 2762 } 2763 ierr = TSViewFromOptions(ts,NULL,"-ts_view");CHKERRQ(ierr); 2764 ierr = PetscObjectSAWsBlock((PetscObject)ts);CHKERRQ(ierr); 2765 PetscFunctionReturn(0); 2766 } 2767 2768 #undef __FUNCT__ 2769 #define __FUNCT__ "TSMonitor" 2770 /*@ 2771 TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet() 2772 2773 Collective on TS 2774 2775 Input Parameters: 2776 + ts - time stepping context obtained from TSCreate() 2777 . step - step number that has just completed 2778 . ptime - model time of the state 2779 - u - state at the current model time 2780 2781 Notes: 2782 TSMonitor() is typically used within the time stepping implementations. 2783 Users might call this function when using the TSStep() interface instead of TSSolve(). 2784 2785 Level: advanced 2786 2787 .keywords: TS, timestep 2788 @*/ 2789 PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u) 2790 { 2791 PetscErrorCode ierr; 2792 PetscInt i,n = ts->numbermonitors; 2793 2794 PetscFunctionBegin; 2795 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2796 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 2797 for (i=0; i<n; i++) { 2798 ierr = (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);CHKERRQ(ierr); 2799 } 2800 PetscFunctionReturn(0); 2801 } 2802 2803 /* ------------------------------------------------------------------------*/ 2804 #undef __FUNCT__ 2805 #define __FUNCT__ "TSMonitorLGCtxCreate" 2806 /*@C 2807 TSMonitorLGCtxCreate - Creates a line graph context for use with 2808 TS to monitor the solution process graphically in various ways 2809 2810 Collective on TS 2811 2812 Input Parameters: 2813 + host - the X display to open, or null for the local machine 2814 . label - the title to put in the title bar 2815 . x, y - the screen coordinates of the upper left coordinate of the window 2816 . m, n - the screen width and height in pixels 2817 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time 2818 2819 Output Parameter: 2820 . ctx - the context 2821 2822 Options Database Key: 2823 + -ts_monitor_lg_timestep - automatically sets line graph monitor 2824 . -ts_monitor_lg_solution - 2825 . -ts_monitor_lg_error - 2826 . -ts_monitor_lg_ksp_iterations - 2827 . -ts_monitor_lg_snes_iterations - 2828 - -lg_indicate_data_points <true,false> - indicate the data points (at each time step) on the plot; default is true 2829 2830 Notes: 2831 Use TSMonitorLGCtxDestroy() to destroy. 2832 2833 Level: intermediate 2834 2835 .keywords: TS, monitor, line graph, residual, seealso 2836 2837 .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 2838 2839 @*/ 2840 PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx) 2841 { 2842 PetscDraw win; 2843 PetscErrorCode ierr; 2844 2845 PetscFunctionBegin; 2846 ierr = PetscNew(ctx);CHKERRQ(ierr); 2847 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 2848 ierr = PetscDrawSetFromOptions(win);CHKERRQ(ierr); 2849 ierr = PetscDrawLGCreate(win,1,&(*ctx)->lg);CHKERRQ(ierr); 2850 ierr = PetscLogObjectParent((PetscObject)(*ctx)->lg,(PetscObject)win);CHKERRQ(ierr); 2851 ierr = PetscDrawLGIndicateDataPoints((*ctx)->lg,PETSC_TRUE);CHKERRQ(ierr); 2852 ierr = PetscDrawLGSetFromOptions((*ctx)->lg);CHKERRQ(ierr); 2853 (*ctx)->howoften = howoften; 2854 PetscFunctionReturn(0); 2855 } 2856 2857 #undef __FUNCT__ 2858 #define __FUNCT__ "TSMonitorLGTimeStep" 2859 PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx) 2860 { 2861 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 2862 PetscReal x = ptime,y; 2863 PetscErrorCode ierr; 2864 2865 PetscFunctionBegin; 2866 if (!step) { 2867 PetscDrawAxis axis; 2868 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 2869 ierr = PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time","Time step");CHKERRQ(ierr); 2870 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 2871 ierr = PetscDrawLGIndicateDataPoints(ctx->lg,PETSC_TRUE);CHKERRQ(ierr); 2872 } 2873 ierr = TSGetTimeStep(ts,&y);CHKERRQ(ierr); 2874 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 2875 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 2876 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 2877 } 2878 PetscFunctionReturn(0); 2879 } 2880 2881 #undef __FUNCT__ 2882 #define __FUNCT__ "TSMonitorLGCtxDestroy" 2883 /*@C 2884 TSMonitorLGCtxDestroy - Destroys a line graph context that was created 2885 with TSMonitorLGCtxCreate(). 2886 2887 Collective on TSMonitorLGCtx 2888 2889 Input Parameter: 2890 . ctx - the monitor context 2891 2892 Level: intermediate 2893 2894 .keywords: TS, monitor, line graph, destroy 2895 2896 .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep(); 2897 @*/ 2898 PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx) 2899 { 2900 PetscDraw draw; 2901 PetscErrorCode ierr; 2902 2903 PetscFunctionBegin; 2904 ierr = PetscDrawLGGetDraw((*ctx)->lg,&draw);CHKERRQ(ierr); 2905 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 2906 ierr = PetscDrawLGDestroy(&(*ctx)->lg);CHKERRQ(ierr); 2907 ierr = PetscFree(*ctx);CHKERRQ(ierr); 2908 PetscFunctionReturn(0); 2909 } 2910 2911 #undef __FUNCT__ 2912 #define __FUNCT__ "TSGetTime" 2913 /*@ 2914 TSGetTime - Gets the time of the most recently completed step. 2915 2916 Not Collective 2917 2918 Input Parameter: 2919 . ts - the TS context obtained from TSCreate() 2920 2921 Output Parameter: 2922 . t - the current time 2923 2924 Level: beginner 2925 2926 Note: 2927 When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(), 2928 TSSetPreStage(), TSSetPostStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated. 2929 2930 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 2931 2932 .keywords: TS, get, time 2933 @*/ 2934 PetscErrorCode TSGetTime(TS ts,PetscReal *t) 2935 { 2936 PetscFunctionBegin; 2937 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2938 PetscValidRealPointer(t,2); 2939 *t = ts->ptime; 2940 PetscFunctionReturn(0); 2941 } 2942 2943 #undef __FUNCT__ 2944 #define __FUNCT__ "TSSetTime" 2945 /*@ 2946 TSSetTime - Allows one to reset the time. 2947 2948 Logically Collective on TS 2949 2950 Input Parameters: 2951 + ts - the TS context obtained from TSCreate() 2952 - time - the time 2953 2954 Level: intermediate 2955 2956 .seealso: TSGetTime(), TSSetDuration() 2957 2958 .keywords: TS, set, time 2959 @*/ 2960 PetscErrorCode TSSetTime(TS ts, PetscReal t) 2961 { 2962 PetscFunctionBegin; 2963 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2964 PetscValidLogicalCollectiveReal(ts,t,2); 2965 ts->ptime = t; 2966 PetscFunctionReturn(0); 2967 } 2968 2969 #undef __FUNCT__ 2970 #define __FUNCT__ "TSSetOptionsPrefix" 2971 /*@C 2972 TSSetOptionsPrefix - Sets the prefix used for searching for all 2973 TS options in the database. 2974 2975 Logically Collective on TS 2976 2977 Input Parameter: 2978 + ts - The TS context 2979 - prefix - The prefix to prepend to all option names 2980 2981 Notes: 2982 A hyphen (-) must NOT be given at the beginning of the prefix name. 2983 The first character of all runtime options is AUTOMATICALLY the 2984 hyphen. 2985 2986 Level: advanced 2987 2988 .keywords: TS, set, options, prefix, database 2989 2990 .seealso: TSSetFromOptions() 2991 2992 @*/ 2993 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 2994 { 2995 PetscErrorCode ierr; 2996 SNES snes; 2997 2998 PetscFunctionBegin; 2999 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3000 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3001 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3002 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 3003 PetscFunctionReturn(0); 3004 } 3005 3006 3007 #undef __FUNCT__ 3008 #define __FUNCT__ "TSAppendOptionsPrefix" 3009 /*@C 3010 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 3011 TS options in the database. 3012 3013 Logically Collective on TS 3014 3015 Input Parameter: 3016 + ts - The TS context 3017 - prefix - The prefix to prepend to all option names 3018 3019 Notes: 3020 A hyphen (-) must NOT be given at the beginning of the prefix name. 3021 The first character of all runtime options is AUTOMATICALLY the 3022 hyphen. 3023 3024 Level: advanced 3025 3026 .keywords: TS, append, options, prefix, database 3027 3028 .seealso: TSGetOptionsPrefix() 3029 3030 @*/ 3031 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 3032 { 3033 PetscErrorCode ierr; 3034 SNES snes; 3035 3036 PetscFunctionBegin; 3037 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3038 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3039 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3040 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 3041 PetscFunctionReturn(0); 3042 } 3043 3044 #undef __FUNCT__ 3045 #define __FUNCT__ "TSGetOptionsPrefix" 3046 /*@C 3047 TSGetOptionsPrefix - Sets the prefix used for searching for all 3048 TS options in the database. 3049 3050 Not Collective 3051 3052 Input Parameter: 3053 . ts - The TS context 3054 3055 Output Parameter: 3056 . prefix - A pointer to the prefix string used 3057 3058 Notes: On the fortran side, the user should pass in a string 'prifix' of 3059 sufficient length to hold the prefix. 3060 3061 Level: intermediate 3062 3063 .keywords: TS, get, options, prefix, database 3064 3065 .seealso: TSAppendOptionsPrefix() 3066 @*/ 3067 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 3068 { 3069 PetscErrorCode ierr; 3070 3071 PetscFunctionBegin; 3072 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3073 PetscValidPointer(prefix,2); 3074 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3075 PetscFunctionReturn(0); 3076 } 3077 3078 #undef __FUNCT__ 3079 #define __FUNCT__ "TSGetRHSJacobian" 3080 /*@C 3081 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 3082 3083 Not Collective, but parallel objects are returned if TS is parallel 3084 3085 Input Parameter: 3086 . ts - The TS context obtained from TSCreate() 3087 3088 Output Parameters: 3089 + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or NULL) 3090 . Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat (or NULL) 3091 . func - Function to compute the Jacobian of the RHS (or NULL) 3092 - ctx - User-defined context for Jacobian evaluation routine (or NULL) 3093 3094 Notes: You can pass in NULL for any return argument you do not need. 3095 3096 Level: intermediate 3097 3098 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3099 3100 .keywords: TS, timestep, get, matrix, Jacobian 3101 @*/ 3102 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx) 3103 { 3104 PetscErrorCode ierr; 3105 SNES snes; 3106 DM dm; 3107 3108 PetscFunctionBegin; 3109 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3110 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3111 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3112 ierr = DMTSGetRHSJacobian(dm,func,ctx);CHKERRQ(ierr); 3113 PetscFunctionReturn(0); 3114 } 3115 3116 #undef __FUNCT__ 3117 #define __FUNCT__ "TSGetIJacobian" 3118 /*@C 3119 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 3120 3121 Not Collective, but parallel objects are returned if TS is parallel 3122 3123 Input Parameter: 3124 . ts - The TS context obtained from TSCreate() 3125 3126 Output Parameters: 3127 + Amat - The (approximate) Jacobian of F(t,U,U_t) 3128 . Pmat - The matrix from which the preconditioner is constructed, often the same as Amat 3129 . f - The function to compute the matrices 3130 - ctx - User-defined context for Jacobian evaluation routine 3131 3132 Notes: You can pass in NULL for any return argument you do not need. 3133 3134 Level: advanced 3135 3136 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3137 3138 .keywords: TS, timestep, get, matrix, Jacobian 3139 @*/ 3140 PetscErrorCode TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx) 3141 { 3142 PetscErrorCode ierr; 3143 SNES snes; 3144 DM dm; 3145 3146 PetscFunctionBegin; 3147 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3148 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 3149 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3150 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3151 ierr = DMTSGetIJacobian(dm,f,ctx);CHKERRQ(ierr); 3152 PetscFunctionReturn(0); 3153 } 3154 3155 3156 #undef __FUNCT__ 3157 #define __FUNCT__ "TSMonitorDrawSolution" 3158 /*@C 3159 TSMonitorDrawSolution - Monitors progress of the TS solvers by calling 3160 VecView() for the solution at each timestep 3161 3162 Collective on TS 3163 3164 Input Parameters: 3165 + ts - the TS context 3166 . step - current time-step 3167 . ptime - current time 3168 - dummy - either a viewer or NULL 3169 3170 Options Database: 3171 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3172 3173 Notes: the initial solution and current solution are not displayed with a common axis scaling so generally the option -ts_monitor_draw_solution_initial 3174 will look bad 3175 3176 Level: intermediate 3177 3178 .keywords: TS, vector, monitor, view 3179 3180 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3181 @*/ 3182 PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3183 { 3184 PetscErrorCode ierr; 3185 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3186 PetscDraw draw; 3187 3188 PetscFunctionBegin; 3189 if (!step && ictx->showinitial) { 3190 if (!ictx->initialsolution) { 3191 ierr = VecDuplicate(u,&ictx->initialsolution);CHKERRQ(ierr); 3192 } 3193 ierr = VecCopy(u,ictx->initialsolution);CHKERRQ(ierr); 3194 } 3195 if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 3196 3197 if (ictx->showinitial) { 3198 PetscReal pause; 3199 ierr = PetscViewerDrawGetPause(ictx->viewer,&pause);CHKERRQ(ierr); 3200 ierr = PetscViewerDrawSetPause(ictx->viewer,0.0);CHKERRQ(ierr); 3201 ierr = VecView(ictx->initialsolution,ictx->viewer);CHKERRQ(ierr); 3202 ierr = PetscViewerDrawSetPause(ictx->viewer,pause);CHKERRQ(ierr); 3203 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);CHKERRQ(ierr); 3204 } 3205 ierr = VecView(u,ictx->viewer);CHKERRQ(ierr); 3206 if (ictx->showtimestepandtime) { 3207 PetscReal xl,yl,xr,yr,tw,w,h; 3208 char time[32]; 3209 size_t len; 3210 3211 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3212 ierr = PetscSNPrintf(time,32,"Timestep %d Time %f",(int)step,(double)ptime);CHKERRQ(ierr); 3213 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3214 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3215 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3216 w = xl + .5*(xr - xl) - .5*len*tw; 3217 h = yl + .95*(yr - yl); 3218 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3219 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3220 } 3221 3222 if (ictx->showinitial) { 3223 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);CHKERRQ(ierr); 3224 } 3225 PetscFunctionReturn(0); 3226 } 3227 3228 #undef __FUNCT__ 3229 #define __FUNCT__ "TSMonitorDrawSolutionPhase" 3230 /*@C 3231 TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram 3232 3233 Collective on TS 3234 3235 Input Parameters: 3236 + ts - the TS context 3237 . step - current time-step 3238 . ptime - current time 3239 - dummy - either a viewer or NULL 3240 3241 Level: intermediate 3242 3243 .keywords: TS, vector, monitor, view 3244 3245 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3246 @*/ 3247 PetscErrorCode TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3248 { 3249 PetscErrorCode ierr; 3250 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3251 PetscDraw draw; 3252 MPI_Comm comm; 3253 PetscInt n; 3254 PetscMPIInt size; 3255 PetscReal xl,yl,xr,yr,tw,w,h; 3256 char time[32]; 3257 size_t len; 3258 const PetscScalar *U; 3259 3260 PetscFunctionBegin; 3261 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 3262 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3263 if (size != 1) SETERRQ(comm,PETSC_ERR_SUP,"Only allowed for sequential runs"); 3264 ierr = VecGetSize(u,&n);CHKERRQ(ierr); 3265 if (n != 2) SETERRQ(comm,PETSC_ERR_SUP,"Only for ODEs with two unknowns"); 3266 3267 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3268 3269 ierr = VecGetArrayRead(u,&U);CHKERRQ(ierr); 3270 ierr = PetscDrawAxisGetLimits(ictx->axis,&xl,&xr,&yl,&yr);CHKERRQ(ierr); 3271 if ((PetscRealPart(U[0]) < xl) || (PetscRealPart(U[1]) < yl) || (PetscRealPart(U[0]) > xr) || (PetscRealPart(U[1]) > yr)) { 3272 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3273 PetscFunctionReturn(0); 3274 } 3275 if (!step) ictx->color++; 3276 ierr = PetscDrawPoint(draw,PetscRealPart(U[0]),PetscRealPart(U[1]),ictx->color);CHKERRQ(ierr); 3277 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3278 3279 if (ictx->showtimestepandtime) { 3280 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3281 ierr = PetscSNPrintf(time,32,"Timestep %d Time %f",(int)step,(double)ptime);CHKERRQ(ierr); 3282 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3283 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3284 w = xl + .5*(xr - xl) - .5*len*tw; 3285 h = yl + .95*(yr - yl); 3286 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3287 } 3288 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3289 PetscFunctionReturn(0); 3290 } 3291 3292 3293 #undef __FUNCT__ 3294 #define __FUNCT__ "TSMonitorDrawCtxDestroy" 3295 /*@C 3296 TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution() 3297 3298 Collective on TS 3299 3300 Input Parameters: 3301 . ctx - the monitor context 3302 3303 Level: intermediate 3304 3305 .keywords: TS, vector, monitor, view 3306 3307 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError() 3308 @*/ 3309 PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx) 3310 { 3311 PetscErrorCode ierr; 3312 3313 PetscFunctionBegin; 3314 ierr = PetscDrawAxisDestroy(&(*ictx)->axis);CHKERRQ(ierr); 3315 ierr = PetscViewerDestroy(&(*ictx)->viewer);CHKERRQ(ierr); 3316 ierr = VecDestroy(&(*ictx)->initialsolution);CHKERRQ(ierr); 3317 ierr = PetscFree(*ictx);CHKERRQ(ierr); 3318 PetscFunctionReturn(0); 3319 } 3320 3321 #undef __FUNCT__ 3322 #define __FUNCT__ "TSMonitorDrawCtxCreate" 3323 /*@C 3324 TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx 3325 3326 Collective on TS 3327 3328 Input Parameter: 3329 . ts - time-step context 3330 3331 Output Patameter: 3332 . ctx - the monitor context 3333 3334 Options Database: 3335 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3336 3337 Level: intermediate 3338 3339 .keywords: TS, vector, monitor, view 3340 3341 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx() 3342 @*/ 3343 PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx) 3344 { 3345 PetscErrorCode ierr; 3346 3347 PetscFunctionBegin; 3348 ierr = PetscNew(ctx);CHKERRQ(ierr); 3349 ierr = PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);CHKERRQ(ierr); 3350 ierr = PetscViewerSetFromOptions((*ctx)->viewer);CHKERRQ(ierr); 3351 3352 (*ctx)->howoften = howoften; 3353 (*ctx)->showinitial = PETSC_FALSE; 3354 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);CHKERRQ(ierr); 3355 3356 (*ctx)->showtimestepandtime = PETSC_FALSE; 3357 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);CHKERRQ(ierr); 3358 (*ctx)->color = PETSC_DRAW_WHITE; 3359 PetscFunctionReturn(0); 3360 } 3361 3362 #undef __FUNCT__ 3363 #define __FUNCT__ "TSMonitorDrawError" 3364 /*@C 3365 TSMonitorDrawError - Monitors progress of the TS solvers by calling 3366 VecView() for the error at each timestep 3367 3368 Collective on TS 3369 3370 Input Parameters: 3371 + ts - the TS context 3372 . step - current time-step 3373 . ptime - current time 3374 - dummy - either a viewer or NULL 3375 3376 Level: intermediate 3377 3378 .keywords: TS, vector, monitor, view 3379 3380 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3381 @*/ 3382 PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3383 { 3384 PetscErrorCode ierr; 3385 TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy; 3386 PetscViewer viewer = ctx->viewer; 3387 Vec work; 3388 3389 PetscFunctionBegin; 3390 if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 3391 ierr = VecDuplicate(u,&work);CHKERRQ(ierr); 3392 ierr = TSComputeSolutionFunction(ts,ptime,work);CHKERRQ(ierr); 3393 ierr = VecAXPY(work,-1.0,u);CHKERRQ(ierr); 3394 ierr = VecView(work,viewer);CHKERRQ(ierr); 3395 ierr = VecDestroy(&work);CHKERRQ(ierr); 3396 PetscFunctionReturn(0); 3397 } 3398 3399 #include <petsc-private/dmimpl.h> 3400 #undef __FUNCT__ 3401 #define __FUNCT__ "TSSetDM" 3402 /*@ 3403 TSSetDM - Sets the DM that may be used by some preconditioners 3404 3405 Logically Collective on TS and DM 3406 3407 Input Parameters: 3408 + ts - the preconditioner context 3409 - dm - the dm 3410 3411 Level: intermediate 3412 3413 3414 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 3415 @*/ 3416 PetscErrorCode TSSetDM(TS ts,DM dm) 3417 { 3418 PetscErrorCode ierr; 3419 SNES snes; 3420 DMTS tsdm; 3421 3422 PetscFunctionBegin; 3423 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3424 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 3425 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 3426 if (ts->dm->dmts && !dm->dmts) { 3427 ierr = DMCopyDMTS(ts->dm,dm);CHKERRQ(ierr); 3428 ierr = DMGetDMTS(ts->dm,&tsdm);CHKERRQ(ierr); 3429 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 3430 tsdm->originaldm = dm; 3431 } 3432 } 3433 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 3434 } 3435 ts->dm = dm; 3436 3437 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3438 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 3439 PetscFunctionReturn(0); 3440 } 3441 3442 #undef __FUNCT__ 3443 #define __FUNCT__ "TSGetDM" 3444 /*@ 3445 TSGetDM - Gets the DM that may be used by some preconditioners 3446 3447 Not Collective 3448 3449 Input Parameter: 3450 . ts - the preconditioner context 3451 3452 Output Parameter: 3453 . dm - the dm 3454 3455 Level: intermediate 3456 3457 3458 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 3459 @*/ 3460 PetscErrorCode TSGetDM(TS ts,DM *dm) 3461 { 3462 PetscErrorCode ierr; 3463 3464 PetscFunctionBegin; 3465 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3466 if (!ts->dm) { 3467 ierr = DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);CHKERRQ(ierr); 3468 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 3469 } 3470 *dm = ts->dm; 3471 PetscFunctionReturn(0); 3472 } 3473 3474 #undef __FUNCT__ 3475 #define __FUNCT__ "SNESTSFormFunction" 3476 /*@ 3477 SNESTSFormFunction - Function to evaluate nonlinear residual 3478 3479 Logically Collective on SNES 3480 3481 Input Parameter: 3482 + snes - nonlinear solver 3483 . U - the current state at which to evaluate the residual 3484 - ctx - user context, must be a TS 3485 3486 Output Parameter: 3487 . F - the nonlinear residual 3488 3489 Notes: 3490 This function is not normally called by users and is automatically registered with the SNES used by TS. 3491 It is most frequently passed to MatFDColoringSetFunction(). 3492 3493 Level: advanced 3494 3495 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 3496 @*/ 3497 PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx) 3498 { 3499 TS ts = (TS)ctx; 3500 PetscErrorCode ierr; 3501 3502 PetscFunctionBegin; 3503 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3504 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3505 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 3506 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 3507 ierr = (ts->ops->snesfunction)(snes,U,F,ts);CHKERRQ(ierr); 3508 PetscFunctionReturn(0); 3509 } 3510 3511 #undef __FUNCT__ 3512 #define __FUNCT__ "SNESTSFormJacobian" 3513 /*@ 3514 SNESTSFormJacobian - Function to evaluate the Jacobian 3515 3516 Collective on SNES 3517 3518 Input Parameter: 3519 + snes - nonlinear solver 3520 . U - the current state at which to evaluate the residual 3521 - ctx - user context, must be a TS 3522 3523 Output Parameter: 3524 + A - the Jacobian 3525 . B - the preconditioning matrix (may be the same as A) 3526 - flag - indicates any structure change in the matrix 3527 3528 Notes: 3529 This function is not normally called by users and is automatically registered with the SNES used by TS. 3530 3531 Level: developer 3532 3533 .seealso: SNESSetJacobian() 3534 @*/ 3535 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat A,Mat B,void *ctx) 3536 { 3537 TS ts = (TS)ctx; 3538 PetscErrorCode ierr; 3539 3540 PetscFunctionBegin; 3541 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3542 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3543 PetscValidPointer(A,3); 3544 PetscValidHeaderSpecific(A,MAT_CLASSID,3); 3545 PetscValidPointer(B,4); 3546 PetscValidHeaderSpecific(B,MAT_CLASSID,4); 3547 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 3548 ierr = (ts->ops->snesjacobian)(snes,U,A,B,ts);CHKERRQ(ierr); 3549 PetscFunctionReturn(0); 3550 } 3551 3552 #undef __FUNCT__ 3553 #define __FUNCT__ "TSComputeRHSFunctionLinear" 3554 /*@C 3555 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 3556 3557 Collective on TS 3558 3559 Input Arguments: 3560 + ts - time stepping context 3561 . t - time at which to evaluate 3562 . U - state at which to evaluate 3563 - ctx - context 3564 3565 Output Arguments: 3566 . F - right hand side 3567 3568 Level: intermediate 3569 3570 Notes: 3571 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 3572 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 3573 3574 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 3575 @*/ 3576 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 3577 { 3578 PetscErrorCode ierr; 3579 Mat Arhs,Brhs; 3580 3581 PetscFunctionBegin; 3582 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 3583 ierr = TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);CHKERRQ(ierr); 3584 ierr = MatMult(Arhs,U,F);CHKERRQ(ierr); 3585 PetscFunctionReturn(0); 3586 } 3587 3588 #undef __FUNCT__ 3589 #define __FUNCT__ "TSComputeRHSJacobianConstant" 3590 /*@C 3591 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 3592 3593 Collective on TS 3594 3595 Input Arguments: 3596 + ts - time stepping context 3597 . t - time at which to evaluate 3598 . U - state at which to evaluate 3599 - ctx - context 3600 3601 Output Arguments: 3602 + A - pointer to operator 3603 . B - pointer to preconditioning matrix 3604 - flg - matrix structure flag 3605 3606 Level: intermediate 3607 3608 Notes: 3609 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 3610 3611 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 3612 @*/ 3613 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx) 3614 { 3615 PetscFunctionBegin; 3616 PetscFunctionReturn(0); 3617 } 3618 3619 #undef __FUNCT__ 3620 #define __FUNCT__ "TSComputeIFunctionLinear" 3621 /*@C 3622 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 3623 3624 Collective on TS 3625 3626 Input Arguments: 3627 + ts - time stepping context 3628 . t - time at which to evaluate 3629 . U - state at which to evaluate 3630 . Udot - time derivative of state vector 3631 - ctx - context 3632 3633 Output Arguments: 3634 . F - left hand side 3635 3636 Level: intermediate 3637 3638 Notes: 3639 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 3640 user is required to write their own TSComputeIFunction. 3641 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 3642 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 3643 3644 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 3645 @*/ 3646 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 3647 { 3648 PetscErrorCode ierr; 3649 Mat A,B; 3650 3651 PetscFunctionBegin; 3652 ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr); 3653 ierr = TSComputeIJacobian(ts,t,U,Udot,1.0,A,B,PETSC_TRUE);CHKERRQ(ierr); 3654 ierr = MatMult(A,Udot,F);CHKERRQ(ierr); 3655 PetscFunctionReturn(0); 3656 } 3657 3658 #undef __FUNCT__ 3659 #define __FUNCT__ "TSComputeIJacobianConstant" 3660 /*@C 3661 TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE 3662 3663 Collective on TS 3664 3665 Input Arguments: 3666 + ts - time stepping context 3667 . t - time at which to evaluate 3668 . U - state at which to evaluate 3669 . Udot - time derivative of state vector 3670 . shift - shift to apply 3671 - ctx - context 3672 3673 Output Arguments: 3674 + A - pointer to operator 3675 . B - pointer to preconditioning matrix 3676 - flg - matrix structure flag 3677 3678 Level: advanced 3679 3680 Notes: 3681 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 3682 3683 It is only appropriate for problems of the form 3684 3685 $ M Udot = F(U,t) 3686 3687 where M is constant and F is non-stiff. The user must pass M to TSSetIJacobian(). The current implementation only 3688 works with IMEX time integration methods such as TSROSW and TSARKIMEX, since there is no support for de-constructing 3689 an implicit operator of the form 3690 3691 $ shift*M + J 3692 3693 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 3694 a copy of M or reassemble it when requested. 3695 3696 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 3697 @*/ 3698 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,void *ctx) 3699 { 3700 PetscErrorCode ierr; 3701 3702 PetscFunctionBegin; 3703 ierr = MatScale(A, shift / ts->ijacobian.shift);CHKERRQ(ierr); 3704 ts->ijacobian.shift = shift; 3705 PetscFunctionReturn(0); 3706 } 3707 3708 #undef __FUNCT__ 3709 #define __FUNCT__ "TSGetEquationType" 3710 /*@ 3711 TSGetEquationType - Gets the type of the equation that TS is solving. 3712 3713 Not Collective 3714 3715 Input Parameter: 3716 . ts - the TS context 3717 3718 Output Parameter: 3719 . equation_type - see TSEquationType 3720 3721 Level: beginner 3722 3723 .keywords: TS, equation type 3724 3725 .seealso: TSSetEquationType(), TSEquationType 3726 @*/ 3727 PetscErrorCode TSGetEquationType(TS ts,TSEquationType *equation_type) 3728 { 3729 PetscFunctionBegin; 3730 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3731 PetscValidPointer(equation_type,2); 3732 *equation_type = ts->equation_type; 3733 PetscFunctionReturn(0); 3734 } 3735 3736 #undef __FUNCT__ 3737 #define __FUNCT__ "TSSetEquationType" 3738 /*@ 3739 TSSetEquationType - Sets the type of the equation that TS is solving. 3740 3741 Not Collective 3742 3743 Input Parameter: 3744 + ts - the TS context 3745 . equation_type - see TSEquationType 3746 3747 Level: advanced 3748 3749 .keywords: TS, equation type 3750 3751 .seealso: TSGetEquationType(), TSEquationType 3752 @*/ 3753 PetscErrorCode TSSetEquationType(TS ts,TSEquationType equation_type) 3754 { 3755 PetscFunctionBegin; 3756 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3757 ts->equation_type = equation_type; 3758 PetscFunctionReturn(0); 3759 } 3760 3761 #undef __FUNCT__ 3762 #define __FUNCT__ "TSGetConvergedReason" 3763 /*@ 3764 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 3765 3766 Not Collective 3767 3768 Input Parameter: 3769 . ts - the TS context 3770 3771 Output Parameter: 3772 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3773 manual pages for the individual convergence tests for complete lists 3774 3775 Level: beginner 3776 3777 Notes: 3778 Can only be called after the call to TSSolve() is complete. 3779 3780 .keywords: TS, nonlinear, set, convergence, test 3781 3782 .seealso: TSSetConvergenceTest(), TSConvergedReason 3783 @*/ 3784 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 3785 { 3786 PetscFunctionBegin; 3787 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3788 PetscValidPointer(reason,2); 3789 *reason = ts->reason; 3790 PetscFunctionReturn(0); 3791 } 3792 3793 #undef __FUNCT__ 3794 #define __FUNCT__ "TSSetConvergedReason" 3795 /*@ 3796 TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve. 3797 3798 Not Collective 3799 3800 Input Parameter: 3801 + ts - the TS context 3802 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3803 manual pages for the individual convergence tests for complete lists 3804 3805 Level: advanced 3806 3807 Notes: 3808 Can only be called during TSSolve() is active. 3809 3810 .keywords: TS, nonlinear, set, convergence, test 3811 3812 .seealso: TSConvergedReason 3813 @*/ 3814 PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason) 3815 { 3816 PetscFunctionBegin; 3817 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3818 ts->reason = reason; 3819 PetscFunctionReturn(0); 3820 } 3821 3822 #undef __FUNCT__ 3823 #define __FUNCT__ "TSGetSolveTime" 3824 /*@ 3825 TSGetSolveTime - Gets the time after a call to TSSolve() 3826 3827 Not Collective 3828 3829 Input Parameter: 3830 . ts - the TS context 3831 3832 Output Parameter: 3833 . ftime - the final time. This time should correspond to the final time set with TSSetDuration() 3834 3835 Level: beginner 3836 3837 Notes: 3838 Can only be called after the call to TSSolve() is complete. 3839 3840 .keywords: TS, nonlinear, set, convergence, test 3841 3842 .seealso: TSSetConvergenceTest(), TSConvergedReason 3843 @*/ 3844 PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime) 3845 { 3846 PetscFunctionBegin; 3847 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3848 PetscValidPointer(ftime,2); 3849 *ftime = ts->solvetime; 3850 PetscFunctionReturn(0); 3851 } 3852 3853 #undef __FUNCT__ 3854 #define __FUNCT__ "TSGetSNESIterations" 3855 /*@ 3856 TSGetSNESIterations - Gets the total number of nonlinear iterations 3857 used by the time integrator. 3858 3859 Not Collective 3860 3861 Input Parameter: 3862 . ts - TS context 3863 3864 Output Parameter: 3865 . nits - number of nonlinear iterations 3866 3867 Notes: 3868 This counter is reset to zero for each successive call to TSSolve(). 3869 3870 Level: intermediate 3871 3872 .keywords: TS, get, number, nonlinear, iterations 3873 3874 .seealso: TSGetKSPIterations() 3875 @*/ 3876 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 3877 { 3878 PetscFunctionBegin; 3879 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3880 PetscValidIntPointer(nits,2); 3881 *nits = ts->snes_its; 3882 PetscFunctionReturn(0); 3883 } 3884 3885 #undef __FUNCT__ 3886 #define __FUNCT__ "TSGetKSPIterations" 3887 /*@ 3888 TSGetKSPIterations - Gets the total number of linear iterations 3889 used by the time integrator. 3890 3891 Not Collective 3892 3893 Input Parameter: 3894 . ts - TS context 3895 3896 Output Parameter: 3897 . lits - number of linear iterations 3898 3899 Notes: 3900 This counter is reset to zero for each successive call to TSSolve(). 3901 3902 Level: intermediate 3903 3904 .keywords: TS, get, number, linear, iterations 3905 3906 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 3907 @*/ 3908 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 3909 { 3910 PetscFunctionBegin; 3911 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3912 PetscValidIntPointer(lits,2); 3913 *lits = ts->ksp_its; 3914 PetscFunctionReturn(0); 3915 } 3916 3917 #undef __FUNCT__ 3918 #define __FUNCT__ "TSGetStepRejections" 3919 /*@ 3920 TSGetStepRejections - Gets the total number of rejected steps. 3921 3922 Not Collective 3923 3924 Input Parameter: 3925 . ts - TS context 3926 3927 Output Parameter: 3928 . rejects - number of steps rejected 3929 3930 Notes: 3931 This counter is reset to zero for each successive call to TSSolve(). 3932 3933 Level: intermediate 3934 3935 .keywords: TS, get, number 3936 3937 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 3938 @*/ 3939 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 3940 { 3941 PetscFunctionBegin; 3942 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3943 PetscValidIntPointer(rejects,2); 3944 *rejects = ts->reject; 3945 PetscFunctionReturn(0); 3946 } 3947 3948 #undef __FUNCT__ 3949 #define __FUNCT__ "TSGetSNESFailures" 3950 /*@ 3951 TSGetSNESFailures - Gets the total number of failed SNES solves 3952 3953 Not Collective 3954 3955 Input Parameter: 3956 . ts - TS context 3957 3958 Output Parameter: 3959 . fails - number of failed nonlinear solves 3960 3961 Notes: 3962 This counter is reset to zero for each successive call to TSSolve(). 3963 3964 Level: intermediate 3965 3966 .keywords: TS, get, number 3967 3968 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 3969 @*/ 3970 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 3971 { 3972 PetscFunctionBegin; 3973 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3974 PetscValidIntPointer(fails,2); 3975 *fails = ts->num_snes_failures; 3976 PetscFunctionReturn(0); 3977 } 3978 3979 #undef __FUNCT__ 3980 #define __FUNCT__ "TSSetMaxStepRejections" 3981 /*@ 3982 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 3983 3984 Not Collective 3985 3986 Input Parameter: 3987 + ts - TS context 3988 - rejects - maximum number of rejected steps, pass -1 for unlimited 3989 3990 Notes: 3991 The counter is reset to zero for each step 3992 3993 Options Database Key: 3994 . -ts_max_reject - Maximum number of step rejections before a step fails 3995 3996 Level: intermediate 3997 3998 .keywords: TS, set, maximum, number 3999 4000 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 4001 @*/ 4002 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 4003 { 4004 PetscFunctionBegin; 4005 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4006 ts->max_reject = rejects; 4007 PetscFunctionReturn(0); 4008 } 4009 4010 #undef __FUNCT__ 4011 #define __FUNCT__ "TSSetMaxSNESFailures" 4012 /*@ 4013 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 4014 4015 Not Collective 4016 4017 Input Parameter: 4018 + ts - TS context 4019 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 4020 4021 Notes: 4022 The counter is reset to zero for each successive call to TSSolve(). 4023 4024 Options Database Key: 4025 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 4026 4027 Level: intermediate 4028 4029 .keywords: TS, set, maximum, number 4030 4031 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 4032 @*/ 4033 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 4034 { 4035 PetscFunctionBegin; 4036 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4037 ts->max_snes_failures = fails; 4038 PetscFunctionReturn(0); 4039 } 4040 4041 #undef __FUNCT__ 4042 #define __FUNCT__ "TSSetErrorIfStepFails" 4043 /*@ 4044 TSSetErrorIfStepFails - Error if no step succeeds 4045 4046 Not Collective 4047 4048 Input Parameter: 4049 + ts - TS context 4050 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 4051 4052 Options Database Key: 4053 . -ts_error_if_step_fails - Error if no step succeeds 4054 4055 Level: intermediate 4056 4057 .keywords: TS, set, error 4058 4059 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 4060 @*/ 4061 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 4062 { 4063 PetscFunctionBegin; 4064 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4065 ts->errorifstepfailed = err; 4066 PetscFunctionReturn(0); 4067 } 4068 4069 #undef __FUNCT__ 4070 #define __FUNCT__ "TSMonitorSolutionBinary" 4071 /*@C 4072 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 4073 4074 Collective on TS 4075 4076 Input Parameters: 4077 + ts - the TS context 4078 . step - current time-step 4079 . ptime - current time 4080 . u - current state 4081 - viewer - binary viewer 4082 4083 Level: intermediate 4084 4085 .keywords: TS, vector, monitor, view 4086 4087 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4088 @*/ 4089 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec u,void *viewer) 4090 { 4091 PetscErrorCode ierr; 4092 PetscViewer v = (PetscViewer)viewer; 4093 4094 PetscFunctionBegin; 4095 ierr = VecView(u,v);CHKERRQ(ierr); 4096 PetscFunctionReturn(0); 4097 } 4098 4099 #undef __FUNCT__ 4100 #define __FUNCT__ "TSMonitorSolutionVTK" 4101 /*@C 4102 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 4103 4104 Collective on TS 4105 4106 Input Parameters: 4107 + ts - the TS context 4108 . step - current time-step 4109 . ptime - current time 4110 . u - current state 4111 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4112 4113 Level: intermediate 4114 4115 Notes: 4116 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. 4117 These are named according to the file name template. 4118 4119 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 4120 4121 .keywords: TS, vector, monitor, view 4122 4123 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4124 @*/ 4125 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate) 4126 { 4127 PetscErrorCode ierr; 4128 char filename[PETSC_MAX_PATH_LEN]; 4129 PetscViewer viewer; 4130 4131 PetscFunctionBegin; 4132 ierr = PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);CHKERRQ(ierr); 4133 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 4134 ierr = VecView(u,viewer);CHKERRQ(ierr); 4135 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4136 PetscFunctionReturn(0); 4137 } 4138 4139 #undef __FUNCT__ 4140 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 4141 /*@C 4142 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 4143 4144 Collective on TS 4145 4146 Input Parameters: 4147 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4148 4149 Level: intermediate 4150 4151 Note: 4152 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 4153 4154 .keywords: TS, vector, monitor, view 4155 4156 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 4157 @*/ 4158 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 4159 { 4160 PetscErrorCode ierr; 4161 4162 PetscFunctionBegin; 4163 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 4164 PetscFunctionReturn(0); 4165 } 4166 4167 #undef __FUNCT__ 4168 #define __FUNCT__ "TSGetAdapt" 4169 /*@ 4170 TSGetAdapt - Get the adaptive controller context for the current method 4171 4172 Collective on TS if controller has not been created yet 4173 4174 Input Arguments: 4175 . ts - time stepping context 4176 4177 Output Arguments: 4178 . adapt - adaptive controller 4179 4180 Level: intermediate 4181 4182 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 4183 @*/ 4184 PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt) 4185 { 4186 PetscErrorCode ierr; 4187 4188 PetscFunctionBegin; 4189 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4190 PetscValidPointer(adapt,2); 4191 if (!ts->adapt) { 4192 ierr = TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);CHKERRQ(ierr); 4193 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->adapt);CHKERRQ(ierr); 4194 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 4195 } 4196 *adapt = ts->adapt; 4197 PetscFunctionReturn(0); 4198 } 4199 4200 #undef __FUNCT__ 4201 #define __FUNCT__ "TSSetTolerances" 4202 /*@ 4203 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 4204 4205 Logically Collective 4206 4207 Input Arguments: 4208 + ts - time integration context 4209 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 4210 . vatol - vector of absolute tolerances or NULL, used in preference to atol if present 4211 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 4212 - vrtol - vector of relative tolerances or NULL, used in preference to atol if present 4213 4214 Level: beginner 4215 4216 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 4217 @*/ 4218 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 4219 { 4220 PetscErrorCode ierr; 4221 4222 PetscFunctionBegin; 4223 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 4224 if (vatol) { 4225 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 4226 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 4227 4228 ts->vatol = vatol; 4229 } 4230 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 4231 if (vrtol) { 4232 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 4233 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 4234 4235 ts->vrtol = vrtol; 4236 } 4237 PetscFunctionReturn(0); 4238 } 4239 4240 #undef __FUNCT__ 4241 #define __FUNCT__ "TSGetTolerances" 4242 /*@ 4243 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 4244 4245 Logically Collective 4246 4247 Input Arguments: 4248 . ts - time integration context 4249 4250 Output Arguments: 4251 + atol - scalar absolute tolerances, NULL to ignore 4252 . vatol - vector of absolute tolerances, NULL to ignore 4253 . rtol - scalar relative tolerances, NULL to ignore 4254 - vrtol - vector of relative tolerances, NULL to ignore 4255 4256 Level: beginner 4257 4258 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 4259 @*/ 4260 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 4261 { 4262 PetscFunctionBegin; 4263 if (atol) *atol = ts->atol; 4264 if (vatol) *vatol = ts->vatol; 4265 if (rtol) *rtol = ts->rtol; 4266 if (vrtol) *vrtol = ts->vrtol; 4267 PetscFunctionReturn(0); 4268 } 4269 4270 #undef __FUNCT__ 4271 #define __FUNCT__ "TSErrorNormWRMS" 4272 /*@ 4273 TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state 4274 4275 Collective on TS 4276 4277 Input Arguments: 4278 + ts - time stepping context 4279 - Y - state vector to be compared to ts->vec_sol 4280 4281 Output Arguments: 4282 . norm - weighted norm, a value of 1.0 is considered small 4283 4284 Level: developer 4285 4286 .seealso: TSSetTolerances() 4287 @*/ 4288 PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm) 4289 { 4290 PetscErrorCode ierr; 4291 PetscInt i,n,N; 4292 const PetscScalar *u,*y; 4293 Vec U; 4294 PetscReal sum,gsum; 4295 4296 PetscFunctionBegin; 4297 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4298 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 4299 PetscValidPointer(norm,3); 4300 U = ts->vec_sol; 4301 PetscCheckSameTypeAndComm(U,1,Y,2); 4302 if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 4303 4304 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 4305 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 4306 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 4307 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 4308 sum = 0.; 4309 if (ts->vatol && ts->vrtol) { 4310 const PetscScalar *atol,*rtol; 4311 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4312 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4313 for (i=0; i<n; i++) { 4314 PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4315 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4316 } 4317 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4318 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4319 } else if (ts->vatol) { /* vector atol, scalar rtol */ 4320 const PetscScalar *atol; 4321 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4322 for (i=0; i<n; i++) { 4323 PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4324 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4325 } 4326 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4327 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 4328 const PetscScalar *rtol; 4329 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4330 for (i=0; i<n; i++) { 4331 PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4332 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4333 } 4334 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4335 } else { /* scalar atol, scalar rtol */ 4336 for (i=0; i<n; i++) { 4337 PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4338 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4339 } 4340 } 4341 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 4342 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 4343 4344 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4345 *norm = PetscSqrtReal(gsum / N); 4346 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 4347 PetscFunctionReturn(0); 4348 } 4349 4350 #undef __FUNCT__ 4351 #define __FUNCT__ "TSSetCFLTimeLocal" 4352 /*@ 4353 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 4354 4355 Logically Collective on TS 4356 4357 Input Arguments: 4358 + ts - time stepping context 4359 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 4360 4361 Note: 4362 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 4363 4364 Level: intermediate 4365 4366 .seealso: TSGetCFLTime(), TSADAPTCFL 4367 @*/ 4368 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 4369 { 4370 PetscFunctionBegin; 4371 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4372 ts->cfltime_local = cfltime; 4373 ts->cfltime = -1.; 4374 PetscFunctionReturn(0); 4375 } 4376 4377 #undef __FUNCT__ 4378 #define __FUNCT__ "TSGetCFLTime" 4379 /*@ 4380 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 4381 4382 Collective on TS 4383 4384 Input Arguments: 4385 . ts - time stepping context 4386 4387 Output Arguments: 4388 . cfltime - maximum stable time step for forward Euler 4389 4390 Level: advanced 4391 4392 .seealso: TSSetCFLTimeLocal() 4393 @*/ 4394 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 4395 { 4396 PetscErrorCode ierr; 4397 4398 PetscFunctionBegin; 4399 if (ts->cfltime < 0) { 4400 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4401 } 4402 *cfltime = ts->cfltime; 4403 PetscFunctionReturn(0); 4404 } 4405 4406 #undef __FUNCT__ 4407 #define __FUNCT__ "TSVISetVariableBounds" 4408 /*@ 4409 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 4410 4411 Input Parameters: 4412 . ts - the TS context. 4413 . xl - lower bound. 4414 . xu - upper bound. 4415 4416 Notes: 4417 If this routine is not called then the lower and upper bounds are set to 4418 PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp(). 4419 4420 Level: advanced 4421 4422 @*/ 4423 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 4424 { 4425 PetscErrorCode ierr; 4426 SNES snes; 4427 4428 PetscFunctionBegin; 4429 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4430 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 4431 PetscFunctionReturn(0); 4432 } 4433 4434 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4435 #include <mex.h> 4436 4437 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 4438 4439 #undef __FUNCT__ 4440 #define __FUNCT__ "TSComputeFunction_Matlab" 4441 /* 4442 TSComputeFunction_Matlab - Calls the function that has been set with 4443 TSSetFunctionMatlab(). 4444 4445 Collective on TS 4446 4447 Input Parameters: 4448 + snes - the TS context 4449 - u - input vector 4450 4451 Output Parameter: 4452 . y - function vector, as set by TSSetFunction() 4453 4454 Notes: 4455 TSComputeFunction() is typically used within nonlinear solvers 4456 implementations, so most users would not generally call this routine 4457 themselves. 4458 4459 Level: developer 4460 4461 .keywords: TS, nonlinear, compute, function 4462 4463 .seealso: TSSetFunction(), TSGetFunction() 4464 */ 4465 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx) 4466 { 4467 PetscErrorCode ierr; 4468 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4469 int nlhs = 1,nrhs = 7; 4470 mxArray *plhs[1],*prhs[7]; 4471 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 4472 4473 PetscFunctionBegin; 4474 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 4475 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4476 PetscValidHeaderSpecific(udot,VEC_CLASSID,4); 4477 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 4478 PetscCheckSameComm(snes,1,u,3); 4479 PetscCheckSameComm(snes,1,y,5); 4480 4481 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4482 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4483 ierr = PetscMemcpy(&lxdot,&udot,sizeof(udot));CHKERRQ(ierr); 4484 ierr = PetscMemcpy(&ly,&y,sizeof(u));CHKERRQ(ierr); 4485 4486 prhs[0] = mxCreateDoubleScalar((double)ls); 4487 prhs[1] = mxCreateDoubleScalar(time); 4488 prhs[2] = mxCreateDoubleScalar((double)lx); 4489 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4490 prhs[4] = mxCreateDoubleScalar((double)ly); 4491 prhs[5] = mxCreateString(sctx->funcname); 4492 prhs[6] = sctx->ctx; 4493 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 4494 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4495 mxDestroyArray(prhs[0]); 4496 mxDestroyArray(prhs[1]); 4497 mxDestroyArray(prhs[2]); 4498 mxDestroyArray(prhs[3]); 4499 mxDestroyArray(prhs[4]); 4500 mxDestroyArray(prhs[5]); 4501 mxDestroyArray(plhs[0]); 4502 PetscFunctionReturn(0); 4503 } 4504 4505 4506 #undef __FUNCT__ 4507 #define __FUNCT__ "TSSetFunctionMatlab" 4508 /* 4509 TSSetFunctionMatlab - Sets the function evaluation routine and function 4510 vector for use by the TS routines in solving ODEs 4511 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4512 4513 Logically Collective on TS 4514 4515 Input Parameters: 4516 + ts - the TS context 4517 - func - function evaluation routine 4518 4519 Calling sequence of func: 4520 $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx); 4521 4522 Level: beginner 4523 4524 .keywords: TS, nonlinear, set, function 4525 4526 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4527 */ 4528 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 4529 { 4530 PetscErrorCode ierr; 4531 TSMatlabContext *sctx; 4532 4533 PetscFunctionBegin; 4534 /* currently sctx is memory bleed */ 4535 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4536 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4537 /* 4538 This should work, but it doesn't 4539 sctx->ctx = ctx; 4540 mexMakeArrayPersistent(sctx->ctx); 4541 */ 4542 sctx->ctx = mxDuplicateArray(ctx); 4543 4544 ierr = TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4545 PetscFunctionReturn(0); 4546 } 4547 4548 #undef __FUNCT__ 4549 #define __FUNCT__ "TSComputeJacobian_Matlab" 4550 /* 4551 TSComputeJacobian_Matlab - Calls the function that has been set with 4552 TSSetJacobianMatlab(). 4553 4554 Collective on TS 4555 4556 Input Parameters: 4557 + ts - the TS context 4558 . u - input vector 4559 . A, B - the matrices 4560 - ctx - user context 4561 4562 Output Parameter: 4563 . flag - structure of the matrix 4564 4565 Level: developer 4566 4567 .keywords: TS, nonlinear, compute, function 4568 4569 .seealso: TSSetFunction(), TSGetFunction() 4570 @*/ 4571 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4572 { 4573 PetscErrorCode ierr; 4574 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4575 int nlhs = 2,nrhs = 9; 4576 mxArray *plhs[2],*prhs[9]; 4577 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 4578 4579 PetscFunctionBegin; 4580 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4581 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4582 4583 /* call Matlab function in ctx with arguments u and y */ 4584 4585 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4586 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4587 ierr = PetscMemcpy(&lxdot,&udot,sizeof(u));CHKERRQ(ierr); 4588 ierr = PetscMemcpy(&lA,A,sizeof(u));CHKERRQ(ierr); 4589 ierr = PetscMemcpy(&lB,B,sizeof(u));CHKERRQ(ierr); 4590 4591 prhs[0] = mxCreateDoubleScalar((double)ls); 4592 prhs[1] = mxCreateDoubleScalar((double)time); 4593 prhs[2] = mxCreateDoubleScalar((double)lx); 4594 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4595 prhs[4] = mxCreateDoubleScalar((double)shift); 4596 prhs[5] = mxCreateDoubleScalar((double)lA); 4597 prhs[6] = mxCreateDoubleScalar((double)lB); 4598 prhs[7] = mxCreateString(sctx->funcname); 4599 prhs[8] = sctx->ctx; 4600 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 4601 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4602 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4603 mxDestroyArray(prhs[0]); 4604 mxDestroyArray(prhs[1]); 4605 mxDestroyArray(prhs[2]); 4606 mxDestroyArray(prhs[3]); 4607 mxDestroyArray(prhs[4]); 4608 mxDestroyArray(prhs[5]); 4609 mxDestroyArray(prhs[6]); 4610 mxDestroyArray(prhs[7]); 4611 mxDestroyArray(plhs[0]); 4612 mxDestroyArray(plhs[1]); 4613 PetscFunctionReturn(0); 4614 } 4615 4616 4617 #undef __FUNCT__ 4618 #define __FUNCT__ "TSSetJacobianMatlab" 4619 /* 4620 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4621 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 4622 4623 Logically Collective on TS 4624 4625 Input Parameters: 4626 + ts - the TS context 4627 . A,B - Jacobian matrices 4628 . func - function evaluation routine 4629 - ctx - user context 4630 4631 Calling sequence of func: 4632 $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx); 4633 4634 4635 Level: developer 4636 4637 .keywords: TS, nonlinear, set, function 4638 4639 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4640 */ 4641 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 4642 { 4643 PetscErrorCode ierr; 4644 TSMatlabContext *sctx; 4645 4646 PetscFunctionBegin; 4647 /* currently sctx is memory bleed */ 4648 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4649 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4650 /* 4651 This should work, but it doesn't 4652 sctx->ctx = ctx; 4653 mexMakeArrayPersistent(sctx->ctx); 4654 */ 4655 sctx->ctx = mxDuplicateArray(ctx); 4656 4657 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4658 PetscFunctionReturn(0); 4659 } 4660 4661 #undef __FUNCT__ 4662 #define __FUNCT__ "TSMonitor_Matlab" 4663 /* 4664 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 4665 4666 Collective on TS 4667 4668 .seealso: TSSetFunction(), TSGetFunction() 4669 @*/ 4670 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx) 4671 { 4672 PetscErrorCode ierr; 4673 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4674 int nlhs = 1,nrhs = 6; 4675 mxArray *plhs[1],*prhs[6]; 4676 long long int lx = 0,ls = 0; 4677 4678 PetscFunctionBegin; 4679 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4680 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 4681 4682 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4683 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4684 4685 prhs[0] = mxCreateDoubleScalar((double)ls); 4686 prhs[1] = mxCreateDoubleScalar((double)it); 4687 prhs[2] = mxCreateDoubleScalar((double)time); 4688 prhs[3] = mxCreateDoubleScalar((double)lx); 4689 prhs[4] = mxCreateString(sctx->funcname); 4690 prhs[5] = sctx->ctx; 4691 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 4692 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4693 mxDestroyArray(prhs[0]); 4694 mxDestroyArray(prhs[1]); 4695 mxDestroyArray(prhs[2]); 4696 mxDestroyArray(prhs[3]); 4697 mxDestroyArray(prhs[4]); 4698 mxDestroyArray(plhs[0]); 4699 PetscFunctionReturn(0); 4700 } 4701 4702 4703 #undef __FUNCT__ 4704 #define __FUNCT__ "TSMonitorSetMatlab" 4705 /* 4706 TSMonitorSetMatlab - Sets the monitor function from Matlab 4707 4708 Level: developer 4709 4710 .keywords: TS, nonlinear, set, function 4711 4712 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4713 */ 4714 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 4715 { 4716 PetscErrorCode ierr; 4717 TSMatlabContext *sctx; 4718 4719 PetscFunctionBegin; 4720 /* currently sctx is memory bleed */ 4721 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4722 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4723 /* 4724 This should work, but it doesn't 4725 sctx->ctx = ctx; 4726 mexMakeArrayPersistent(sctx->ctx); 4727 */ 4728 sctx->ctx = mxDuplicateArray(ctx); 4729 4730 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 4731 PetscFunctionReturn(0); 4732 } 4733 #endif 4734 4735 4736 4737 #undef __FUNCT__ 4738 #define __FUNCT__ "TSMonitorLGSolution" 4739 /*@C 4740 TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector 4741 in a time based line graph 4742 4743 Collective on TS 4744 4745 Input Parameters: 4746 + ts - the TS context 4747 . step - current time-step 4748 . ptime - current time 4749 - lg - a line graph object 4750 4751 Level: intermediate 4752 4753 Notes: each process in a parallel run displays its component solutions in a separate window 4754 4755 .keywords: TS, vector, monitor, view 4756 4757 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4758 @*/ 4759 PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4760 { 4761 PetscErrorCode ierr; 4762 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4763 const PetscScalar *yy; 4764 PetscInt dim; 4765 4766 PetscFunctionBegin; 4767 if (!step) { 4768 PetscDrawAxis axis; 4769 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4770 ierr = PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");CHKERRQ(ierr); 4771 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4772 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4773 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4774 } 4775 ierr = VecGetArrayRead(u,&yy);CHKERRQ(ierr); 4776 #if defined(PETSC_USE_COMPLEX) 4777 { 4778 PetscReal *yreal; 4779 PetscInt i,n; 4780 ierr = VecGetLocalSize(u,&n);CHKERRQ(ierr); 4781 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 4782 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4783 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4784 ierr = PetscFree(yreal);CHKERRQ(ierr); 4785 } 4786 #else 4787 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4788 #endif 4789 ierr = VecRestoreArrayRead(u,&yy);CHKERRQ(ierr); 4790 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 4791 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4792 } 4793 PetscFunctionReturn(0); 4794 } 4795 4796 #undef __FUNCT__ 4797 #define __FUNCT__ "TSMonitorLGError" 4798 /*@C 4799 TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector 4800 in a time based line graph 4801 4802 Collective on TS 4803 4804 Input Parameters: 4805 + ts - the TS context 4806 . step - current time-step 4807 . ptime - current time 4808 - lg - a line graph object 4809 4810 Level: intermediate 4811 4812 Notes: 4813 Only for sequential solves. 4814 4815 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 4816 4817 Options Database Keys: 4818 . -ts_monitor_lg_error - create a graphical monitor of error history 4819 4820 .keywords: TS, vector, monitor, view 4821 4822 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 4823 @*/ 4824 PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4825 { 4826 PetscErrorCode ierr; 4827 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4828 const PetscScalar *yy; 4829 Vec y; 4830 PetscInt dim; 4831 4832 PetscFunctionBegin; 4833 if (!step) { 4834 PetscDrawAxis axis; 4835 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4836 ierr = PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");CHKERRQ(ierr); 4837 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4838 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4839 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4840 } 4841 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 4842 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 4843 ierr = VecAXPY(y,-1.0,u);CHKERRQ(ierr); 4844 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 4845 #if defined(PETSC_USE_COMPLEX) 4846 { 4847 PetscReal *yreal; 4848 PetscInt i,n; 4849 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 4850 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 4851 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4852 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4853 ierr = PetscFree(yreal);CHKERRQ(ierr); 4854 } 4855 #else 4856 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4857 #endif 4858 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 4859 ierr = VecDestroy(&y);CHKERRQ(ierr); 4860 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 4861 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4862 } 4863 PetscFunctionReturn(0); 4864 } 4865 4866 #undef __FUNCT__ 4867 #define __FUNCT__ "TSMonitorLGSNESIterations" 4868 PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4869 { 4870 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4871 PetscReal x = ptime,y; 4872 PetscErrorCode ierr; 4873 PetscInt its; 4874 4875 PetscFunctionBegin; 4876 if (!n) { 4877 PetscDrawAxis axis; 4878 4879 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4880 ierr = PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");CHKERRQ(ierr); 4881 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4882 4883 ctx->snes_its = 0; 4884 } 4885 ierr = TSGetSNESIterations(ts,&its);CHKERRQ(ierr); 4886 y = its - ctx->snes_its; 4887 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4888 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 4889 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4890 } 4891 ctx->snes_its = its; 4892 PetscFunctionReturn(0); 4893 } 4894 4895 #undef __FUNCT__ 4896 #define __FUNCT__ "TSMonitorLGKSPIterations" 4897 PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4898 { 4899 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4900 PetscReal x = ptime,y; 4901 PetscErrorCode ierr; 4902 PetscInt its; 4903 4904 PetscFunctionBegin; 4905 if (!n) { 4906 PetscDrawAxis axis; 4907 4908 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4909 ierr = PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");CHKERRQ(ierr); 4910 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4911 4912 ctx->ksp_its = 0; 4913 } 4914 ierr = TSGetKSPIterations(ts,&its);CHKERRQ(ierr); 4915 y = its - ctx->ksp_its; 4916 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4917 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 4918 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4919 } 4920 ctx->ksp_its = its; 4921 PetscFunctionReturn(0); 4922 } 4923 4924 #undef __FUNCT__ 4925 #define __FUNCT__ "TSComputeLinearStability" 4926 /*@ 4927 TSComputeLinearStability - computes the linear stability function at a point 4928 4929 Collective on TS and Vec 4930 4931 Input Parameters: 4932 + ts - the TS context 4933 - xr,xi - real and imaginary part of input arguments 4934 4935 Output Parameters: 4936 . yr,yi - real and imaginary part of function value 4937 4938 Level: developer 4939 4940 .keywords: TS, compute 4941 4942 .seealso: TSSetRHSFunction(), TSComputeIFunction() 4943 @*/ 4944 PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 4945 { 4946 PetscErrorCode ierr; 4947 4948 PetscFunctionBegin; 4949 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4950 if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method"); 4951 ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr); 4952 PetscFunctionReturn(0); 4953 } 4954 4955 #undef __FUNCT__ 4956 #define __FUNCT__ "TSRollBack" 4957 /*@ 4958 TSRollBack - Rolls back one time step 4959 4960 Collective on TS 4961 4962 Input Parameter: 4963 . ts - the TS context obtained from TSCreate() 4964 4965 Level: advanced 4966 4967 .keywords: TS, timestep, rollback 4968 4969 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate() 4970 @*/ 4971 PetscErrorCode TSRollBack(TS ts) 4972 { 4973 PetscErrorCode ierr; 4974 4975 PetscFunctionBegin; 4976 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 4977 4978 if (!ts->ops->rollback) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRollBack not implemented for type '%s'",((PetscObject)ts)->type_name); 4979 ierr = (*ts->ops->rollback)(ts);CHKERRQ(ierr); 4980 ts->time_step = ts->ptime - ts->ptime_prev; 4981 ts->ptime = ts->ptime_prev; 4982 PetscFunctionReturn(0); 4983 } 4984 4985