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 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 1913 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 1914 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 1915 1916 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 1917 PetscFunctionReturn(0); 1918 } 1919 1920 #undef __FUNCT__ 1921 #define __FUNCT__ "TSGetSNES" 1922 /*@ 1923 TSGetSNES - Returns the SNES (nonlinear solver) associated with 1924 a TS (timestepper) context. Valid only for nonlinear problems. 1925 1926 Not Collective, but SNES is parallel if TS is parallel 1927 1928 Input Parameter: 1929 . ts - the TS context obtained from TSCreate() 1930 1931 Output Parameter: 1932 . snes - the nonlinear solver context 1933 1934 Notes: 1935 The user can then directly manipulate the SNES context to set various 1936 options, etc. Likewise, the user can then extract and manipulate the 1937 KSP, KSP, and PC contexts as well. 1938 1939 TSGetSNES() does not work for integrators that do not use SNES; in 1940 this case TSGetSNES() returns NULL in snes. 1941 1942 Level: beginner 1943 1944 .keywords: timestep, get, SNES 1945 @*/ 1946 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 1947 { 1948 PetscErrorCode ierr; 1949 1950 PetscFunctionBegin; 1951 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1952 PetscValidPointer(snes,2); 1953 if (!ts->snes) { 1954 ierr = SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);CHKERRQ(ierr); 1955 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 1956 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->snes);CHKERRQ(ierr); 1957 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 1958 if (ts->dm) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 1959 if (ts->problem_type == TS_LINEAR) { 1960 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 1961 } 1962 } 1963 *snes = ts->snes; 1964 PetscFunctionReturn(0); 1965 } 1966 1967 #undef __FUNCT__ 1968 #define __FUNCT__ "TSSetSNES" 1969 /*@ 1970 TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context 1971 1972 Collective 1973 1974 Input Parameter: 1975 + ts - the TS context obtained from TSCreate() 1976 - snes - the nonlinear solver context 1977 1978 Notes: 1979 Most users should have the TS created by calling TSGetSNES() 1980 1981 Level: developer 1982 1983 .keywords: timestep, set, SNES 1984 @*/ 1985 PetscErrorCode TSSetSNES(TS ts,SNES snes) 1986 { 1987 PetscErrorCode ierr; 1988 PetscErrorCode (*func)(SNES,Vec,Mat,Mat,void*); 1989 1990 PetscFunctionBegin; 1991 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1992 PetscValidHeaderSpecific(snes,SNES_CLASSID,2); 1993 ierr = PetscObjectReference((PetscObject)snes);CHKERRQ(ierr); 1994 ierr = SNESDestroy(&ts->snes);CHKERRQ(ierr); 1995 1996 ts->snes = snes; 1997 1998 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 1999 ierr = SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);CHKERRQ(ierr); 2000 if (func == SNESTSFormJacobian) { 2001 ierr = SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);CHKERRQ(ierr); 2002 } 2003 PetscFunctionReturn(0); 2004 } 2005 2006 #undef __FUNCT__ 2007 #define __FUNCT__ "TSGetKSP" 2008 /*@ 2009 TSGetKSP - Returns the KSP (linear solver) associated with 2010 a TS (timestepper) context. 2011 2012 Not Collective, but KSP is parallel if TS is parallel 2013 2014 Input Parameter: 2015 . ts - the TS context obtained from TSCreate() 2016 2017 Output Parameter: 2018 . ksp - the nonlinear solver context 2019 2020 Notes: 2021 The user can then directly manipulate the KSP context to set various 2022 options, etc. Likewise, the user can then extract and manipulate the 2023 KSP and PC contexts as well. 2024 2025 TSGetKSP() does not work for integrators that do not use KSP; 2026 in this case TSGetKSP() returns NULL in ksp. 2027 2028 Level: beginner 2029 2030 .keywords: timestep, get, KSP 2031 @*/ 2032 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 2033 { 2034 PetscErrorCode ierr; 2035 SNES snes; 2036 2037 PetscFunctionBegin; 2038 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2039 PetscValidPointer(ksp,2); 2040 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 2041 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 2042 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2043 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 2044 PetscFunctionReturn(0); 2045 } 2046 2047 /* ----------- Routines to set solver parameters ---------- */ 2048 2049 #undef __FUNCT__ 2050 #define __FUNCT__ "TSGetDuration" 2051 /*@ 2052 TSGetDuration - Gets the maximum number of timesteps to use and 2053 maximum time for iteration. 2054 2055 Not Collective 2056 2057 Input Parameters: 2058 + ts - the TS context obtained from TSCreate() 2059 . maxsteps - maximum number of iterations to use, or NULL 2060 - maxtime - final time to iterate to, or NULL 2061 2062 Level: intermediate 2063 2064 .keywords: TS, timestep, get, maximum, iterations, time 2065 @*/ 2066 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 2067 { 2068 PetscFunctionBegin; 2069 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2070 if (maxsteps) { 2071 PetscValidIntPointer(maxsteps,2); 2072 *maxsteps = ts->max_steps; 2073 } 2074 if (maxtime) { 2075 PetscValidScalarPointer(maxtime,3); 2076 *maxtime = ts->max_time; 2077 } 2078 PetscFunctionReturn(0); 2079 } 2080 2081 #undef __FUNCT__ 2082 #define __FUNCT__ "TSSetDuration" 2083 /*@ 2084 TSSetDuration - Sets the maximum number of timesteps to use and 2085 maximum time for iteration. 2086 2087 Logically Collective on TS 2088 2089 Input Parameters: 2090 + ts - the TS context obtained from TSCreate() 2091 . maxsteps - maximum number of iterations to use 2092 - maxtime - final time to iterate to 2093 2094 Options Database Keys: 2095 . -ts_max_steps <maxsteps> - Sets maxsteps 2096 . -ts_final_time <maxtime> - Sets maxtime 2097 2098 Notes: 2099 The default maximum number of iterations is 5000. Default time is 5.0 2100 2101 Level: intermediate 2102 2103 .keywords: TS, timestep, set, maximum, iterations 2104 2105 .seealso: TSSetExactFinalTime() 2106 @*/ 2107 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 2108 { 2109 PetscFunctionBegin; 2110 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2111 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 2112 PetscValidLogicalCollectiveReal(ts,maxtime,2); 2113 if (maxsteps >= 0) ts->max_steps = maxsteps; 2114 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 2115 PetscFunctionReturn(0); 2116 } 2117 2118 #undef __FUNCT__ 2119 #define __FUNCT__ "TSSetSolution" 2120 /*@ 2121 TSSetSolution - Sets the initial solution vector 2122 for use by the TS routines. 2123 2124 Logically Collective on TS and Vec 2125 2126 Input Parameters: 2127 + ts - the TS context obtained from TSCreate() 2128 - u - the solution vector 2129 2130 Level: beginner 2131 2132 .keywords: TS, timestep, set, solution, initial conditions 2133 @*/ 2134 PetscErrorCode TSSetSolution(TS ts,Vec u) 2135 { 2136 PetscErrorCode ierr; 2137 DM dm; 2138 2139 PetscFunctionBegin; 2140 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2141 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2142 ierr = PetscObjectReference((PetscObject)u);CHKERRQ(ierr); 2143 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 2144 2145 ts->vec_sol = u; 2146 2147 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2148 ierr = DMShellSetGlobalVector(dm,u);CHKERRQ(ierr); 2149 PetscFunctionReturn(0); 2150 } 2151 2152 #undef __FUNCT__ 2153 #define __FUNCT__ "TSSetPreStep" 2154 /*@C 2155 TSSetPreStep - Sets the general-purpose function 2156 called once at the beginning of each time step. 2157 2158 Logically Collective on TS 2159 2160 Input Parameters: 2161 + ts - The TS context obtained from TSCreate() 2162 - func - The function 2163 2164 Calling sequence of func: 2165 . func (TS ts); 2166 2167 Level: intermediate 2168 2169 Note: 2170 If a step is rejected, TSStep() will call this routine again before each attempt. 2171 The last completed time step number can be queried using TSGetTimeStepNumber(), the 2172 size of the step being attempted can be obtained using TSGetTimeStep(). 2173 2174 .keywords: TS, timestep 2175 .seealso: TSSetPreStage(), TSSetPostStage(), TSSetPostStep(), TSStep() 2176 @*/ 2177 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 2178 { 2179 PetscFunctionBegin; 2180 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2181 ts->prestep = func; 2182 PetscFunctionReturn(0); 2183 } 2184 2185 #undef __FUNCT__ 2186 #define __FUNCT__ "TSPreStep" 2187 /*@ 2188 TSPreStep - Runs the user-defined pre-step function. 2189 2190 Collective on TS 2191 2192 Input Parameters: 2193 . ts - The TS context obtained from TSCreate() 2194 2195 Notes: 2196 TSPreStep() is typically used within time stepping implementations, 2197 so most users would not generally call this routine themselves. 2198 2199 Level: developer 2200 2201 .keywords: TS, timestep 2202 .seealso: TSSetPreStep(), TSPreStage(), TSPostStage(), TSPostStep() 2203 @*/ 2204 PetscErrorCode TSPreStep(TS ts) 2205 { 2206 PetscErrorCode ierr; 2207 2208 PetscFunctionBegin; 2209 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2210 if (ts->prestep) { 2211 PetscStackCallStandard((*ts->prestep),(ts)); 2212 } 2213 PetscFunctionReturn(0); 2214 } 2215 2216 #undef __FUNCT__ 2217 #define __FUNCT__ "TSSetPreStage" 2218 /*@C 2219 TSSetPreStage - Sets the general-purpose function 2220 called once at the beginning of each stage. 2221 2222 Logically Collective on TS 2223 2224 Input Parameters: 2225 + ts - The TS context obtained from TSCreate() 2226 - func - The function 2227 2228 Calling sequence of func: 2229 . PetscErrorCode func(TS ts, PetscReal stagetime); 2230 2231 Level: intermediate 2232 2233 Note: 2234 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2235 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2236 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2237 2238 .keywords: TS, timestep 2239 .seealso: TSSetPostStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2240 @*/ 2241 PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal)) 2242 { 2243 PetscFunctionBegin; 2244 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2245 ts->prestage = func; 2246 PetscFunctionReturn(0); 2247 } 2248 2249 #undef __FUNCT__ 2250 #define __FUNCT__ "TSSetPostStage" 2251 /*@C 2252 TSSetPostStage - Sets the general-purpose function 2253 called once at the end of each stage. 2254 2255 Logically Collective on TS 2256 2257 Input Parameters: 2258 + ts - The TS context obtained from TSCreate() 2259 - func - The function 2260 2261 Calling sequence of func: 2262 . PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y); 2263 2264 Level: intermediate 2265 2266 Note: 2267 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2268 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2269 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2270 2271 .keywords: TS, timestep 2272 .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2273 @*/ 2274 PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS,PetscReal,PetscInt,Vec*)) 2275 { 2276 PetscFunctionBegin; 2277 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2278 ts->poststage = func; 2279 PetscFunctionReturn(0); 2280 } 2281 2282 #undef __FUNCT__ 2283 #define __FUNCT__ "TSPreStage" 2284 /*@ 2285 TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage() 2286 2287 Collective on TS 2288 2289 Input Parameters: 2290 . ts - The TS context obtained from TSCreate() 2291 stagetime - The absolute time of the current stage 2292 2293 Notes: 2294 TSPreStage() is typically used within time stepping implementations, 2295 most users would not generally call this routine themselves. 2296 2297 Level: developer 2298 2299 .keywords: TS, timestep 2300 .seealso: TSPostStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 2301 @*/ 2302 PetscErrorCode TSPreStage(TS ts, PetscReal stagetime) 2303 { 2304 PetscErrorCode ierr; 2305 2306 PetscFunctionBegin; 2307 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2308 if (ts->prestage) { 2309 PetscStackCallStandard((*ts->prestage),(ts,stagetime)); 2310 } 2311 PetscFunctionReturn(0); 2312 } 2313 2314 #undef __FUNCT__ 2315 #define __FUNCT__ "TSPostStage" 2316 /*@ 2317 TSPostStage - Runs the user-defined post-stage function set using TSSetPostStage() 2318 2319 Collective on TS 2320 2321 Input Parameters: 2322 . ts - The TS context obtained from TSCreate() 2323 stagetime - The absolute time of the current stage 2324 stageindex - Stage number 2325 Y - Array of vectors (of size = total number 2326 of stages) with the stage solutions 2327 2328 Notes: 2329 TSPostStage() is typically used within time stepping implementations, 2330 most users would not generally call this routine themselves. 2331 2332 Level: developer 2333 2334 .keywords: TS, timestep 2335 .seealso: TSPreStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 2336 @*/ 2337 PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y) 2338 { 2339 PetscErrorCode ierr; 2340 2341 PetscFunctionBegin; 2342 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2343 if (ts->prestage) { 2344 PetscStackCallStandard((*ts->poststage),(ts,stagetime,stageindex,Y)); 2345 } 2346 PetscFunctionReturn(0); 2347 } 2348 2349 #undef __FUNCT__ 2350 #define __FUNCT__ "TSSetPostStep" 2351 /*@C 2352 TSSetPostStep - Sets the general-purpose function 2353 called once at the end of each time step. 2354 2355 Logically Collective on TS 2356 2357 Input Parameters: 2358 + ts - The TS context obtained from TSCreate() 2359 - func - The function 2360 2361 Calling sequence of func: 2362 $ func (TS ts); 2363 2364 Level: intermediate 2365 2366 .keywords: TS, timestep 2367 .seealso: TSSetPreStep(), TSSetPreStage(), TSGetTimeStep(), TSGetTimeStepNumber(), TSGetTime() 2368 @*/ 2369 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 2370 { 2371 PetscFunctionBegin; 2372 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2373 ts->poststep = func; 2374 PetscFunctionReturn(0); 2375 } 2376 2377 #undef __FUNCT__ 2378 #define __FUNCT__ "TSPostStep" 2379 /*@ 2380 TSPostStep - Runs the user-defined post-step function. 2381 2382 Collective on TS 2383 2384 Input Parameters: 2385 . ts - The TS context obtained from TSCreate() 2386 2387 Notes: 2388 TSPostStep() is typically used within time stepping implementations, 2389 so most users would not generally call this routine themselves. 2390 2391 Level: developer 2392 2393 .keywords: TS, timestep 2394 @*/ 2395 PetscErrorCode TSPostStep(TS ts) 2396 { 2397 PetscErrorCode ierr; 2398 2399 PetscFunctionBegin; 2400 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2401 if (ts->poststep) { 2402 PetscStackCallStandard((*ts->poststep),(ts)); 2403 } 2404 PetscFunctionReturn(0); 2405 } 2406 2407 /* ------------ Routines to set performance monitoring options ----------- */ 2408 2409 #undef __FUNCT__ 2410 #define __FUNCT__ "TSMonitorSet" 2411 /*@C 2412 TSMonitorSet - Sets an ADDITIONAL function that is to be used at every 2413 timestep to display the iteration's progress. 2414 2415 Logically Collective on TS 2416 2417 Input Parameters: 2418 + ts - the TS context obtained from TSCreate() 2419 . monitor - monitoring routine 2420 . mctx - [optional] user-defined context for private data for the 2421 monitor routine (use NULL if no context is desired) 2422 - monitordestroy - [optional] routine that frees monitor context 2423 (may be NULL) 2424 2425 Calling sequence of monitor: 2426 $ int monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx) 2427 2428 + ts - the TS context 2429 . 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 2430 been interpolated to) 2431 . time - current time 2432 . u - current iterate 2433 - mctx - [optional] monitoring context 2434 2435 Notes: 2436 This routine adds an additional monitor to the list of monitors that 2437 already has been loaded. 2438 2439 Fortran notes: Only a single monitor function can be set for each TS object 2440 2441 Level: intermediate 2442 2443 .keywords: TS, timestep, set, monitor 2444 2445 .seealso: TSMonitorDefault(), TSMonitorCancel() 2446 @*/ 2447 PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**)) 2448 { 2449 PetscFunctionBegin; 2450 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2451 if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 2452 ts->monitor[ts->numbermonitors] = monitor; 2453 ts->monitordestroy[ts->numbermonitors] = mdestroy; 2454 ts->monitorcontext[ts->numbermonitors++] = (void*)mctx; 2455 PetscFunctionReturn(0); 2456 } 2457 2458 #undef __FUNCT__ 2459 #define __FUNCT__ "TSMonitorCancel" 2460 /*@C 2461 TSMonitorCancel - Clears all the monitors that have been set on a time-step object. 2462 2463 Logically Collective on TS 2464 2465 Input Parameters: 2466 . ts - the TS context obtained from TSCreate() 2467 2468 Notes: 2469 There is no way to remove a single, specific monitor. 2470 2471 Level: intermediate 2472 2473 .keywords: TS, timestep, set, monitor 2474 2475 .seealso: TSMonitorDefault(), TSMonitorSet() 2476 @*/ 2477 PetscErrorCode TSMonitorCancel(TS ts) 2478 { 2479 PetscErrorCode ierr; 2480 PetscInt i; 2481 2482 PetscFunctionBegin; 2483 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2484 for (i=0; i<ts->numbermonitors; i++) { 2485 if (ts->monitordestroy[i]) { 2486 ierr = (*ts->monitordestroy[i])(&ts->monitorcontext[i]);CHKERRQ(ierr); 2487 } 2488 } 2489 ts->numbermonitors = 0; 2490 PetscFunctionReturn(0); 2491 } 2492 2493 #undef __FUNCT__ 2494 #define __FUNCT__ "TSMonitorDefault" 2495 /*@ 2496 TSMonitorDefault - Sets the Default monitor 2497 2498 Level: intermediate 2499 2500 .keywords: TS, set, monitor 2501 2502 .seealso: TSMonitorDefault(), TSMonitorSet() 2503 @*/ 2504 PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy) 2505 { 2506 PetscErrorCode ierr; 2507 PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ts)); 2508 2509 PetscFunctionBegin; 2510 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2511 ierr = PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g\n",step,(double)ts->time_step,(double)ptime);CHKERRQ(ierr); 2512 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2513 PetscFunctionReturn(0); 2514 } 2515 2516 #undef __FUNCT__ 2517 #define __FUNCT__ "TSSetRetainStages" 2518 /*@ 2519 TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available. 2520 2521 Logically Collective on TS 2522 2523 Input Argument: 2524 . ts - time stepping context 2525 2526 Output Argument: 2527 . flg - PETSC_TRUE or PETSC_FALSE 2528 2529 Level: intermediate 2530 2531 .keywords: TS, set 2532 2533 .seealso: TSInterpolate(), TSSetPostStep() 2534 @*/ 2535 PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg) 2536 { 2537 PetscFunctionBegin; 2538 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2539 ts->retain_stages = flg; 2540 PetscFunctionReturn(0); 2541 } 2542 2543 #undef __FUNCT__ 2544 #define __FUNCT__ "TSInterpolate" 2545 /*@ 2546 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 2547 2548 Collective on TS 2549 2550 Input Argument: 2551 + ts - time stepping context 2552 - t - time to interpolate to 2553 2554 Output Argument: 2555 . U - state at given time 2556 2557 Notes: 2558 The user should call TSSetRetainStages() before taking a step in which interpolation will be requested. 2559 2560 Level: intermediate 2561 2562 Developer Notes: 2563 TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 2564 2565 .keywords: TS, set 2566 2567 .seealso: TSSetRetainStages(), TSSetPostStep() 2568 @*/ 2569 PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U) 2570 { 2571 PetscErrorCode ierr; 2572 2573 PetscFunctionBegin; 2574 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2575 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 2576 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); 2577 if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name); 2578 ierr = (*ts->ops->interpolate)(ts,t,U);CHKERRQ(ierr); 2579 PetscFunctionReturn(0); 2580 } 2581 2582 #undef __FUNCT__ 2583 #define __FUNCT__ "TSStep" 2584 /*@ 2585 TSStep - Steps one time step 2586 2587 Collective on TS 2588 2589 Input Parameter: 2590 . ts - the TS context obtained from TSCreate() 2591 2592 Level: intermediate 2593 2594 Notes: 2595 The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may 2596 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 2597 2598 This may over-step the final time provided in TSSetDuration() depending on the time-step used. TSSolve() interpolates to exactly the 2599 time provided in TSSetDuration(). One can use TSInterpolate() to determine an interpolated solution within the final timestep. 2600 2601 .keywords: TS, timestep, solve 2602 2603 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSInterpolate() 2604 @*/ 2605 PetscErrorCode TSStep(TS ts) 2606 { 2607 DM dm; 2608 PetscReal ptime_prev; 2609 PetscErrorCode ierr; 2610 static PetscBool cite = PETSC_FALSE; 2611 2612 PetscFunctionBegin; 2613 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2614 ierr = PetscCitationsRegister("@techreport{tspaper,\n" 2615 " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n" 2616 " author = {Shrirang Abhyankar and Jed Brown and Emil Constantinescu and Debojyoti Ghosh and Barry F. Smith},\n" 2617 " type = {Preprint},\n" 2618 " number = {ANL/MCS-P5061-0114},\n" 2619 " institution = {Argonne National Laboratory},\n" 2620 " year = {2014}\n}\n",&cite); 2621 2622 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 2623 ierr = TSSetUp(ts);CHKERRQ(ierr); 2624 2625 ts->reason = TS_CONVERGED_ITERATING; 2626 ptime_prev = ts->ptime; 2627 ierr = DMSetOutputSequenceNumber(dm, ts->steps);CHKERRQ(ierr); 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 - ptime_prev; 2634 ierr = DMSetOutputSequenceNumber(dm, ts->steps);CHKERRQ(ierr); 2635 2636 if (ts->reason < 0) { 2637 if (ts->errorifstepfailed) { 2638 if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) { 2639 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]); 2640 } else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 2641 } 2642 } else if (!ts->reason) { 2643 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 2644 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 2645 } 2646 PetscFunctionReturn(0); 2647 } 2648 2649 #undef __FUNCT__ 2650 #define __FUNCT__ "TSEvaluateStep" 2651 /*@ 2652 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 2653 2654 Collective on TS 2655 2656 Input Arguments: 2657 + ts - time stepping context 2658 . order - desired order of accuracy 2659 - done - whether the step was evaluated at this order (pass NULL to generate an error if not available) 2660 2661 Output Arguments: 2662 . U - state at the end of the current step 2663 2664 Level: advanced 2665 2666 Notes: 2667 This function cannot be called until all stages have been evaluated. 2668 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. 2669 2670 .seealso: TSStep(), TSAdapt 2671 @*/ 2672 PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done) 2673 { 2674 PetscErrorCode ierr; 2675 2676 PetscFunctionBegin; 2677 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2678 PetscValidType(ts,1); 2679 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 2680 if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name); 2681 ierr = (*ts->ops->evaluatestep)(ts,order,U,done);CHKERRQ(ierr); 2682 PetscFunctionReturn(0); 2683 } 2684 2685 #undef __FUNCT__ 2686 #define __FUNCT__ "TSSolve" 2687 /*@ 2688 TSSolve - Steps the requested number of timesteps. 2689 2690 Collective on TS 2691 2692 Input Parameter: 2693 + ts - the TS context obtained from TSCreate() 2694 - u - the solution vector (can be null if TSSetSolution() was used, otherwise must contain the initial conditions) 2695 2696 Level: beginner 2697 2698 Notes: 2699 The final time returned by this function may be different from the time of the internally 2700 held state accessible by TSGetSolution() and TSGetTime() because the method may have 2701 stepped over the final time. 2702 2703 .keywords: TS, timestep, solve 2704 2705 .seealso: TSCreate(), TSSetSolution(), TSStep() 2706 @*/ 2707 PetscErrorCode TSSolve(TS ts,Vec u) 2708 { 2709 Vec solution; 2710 PetscErrorCode ierr; 2711 2712 PetscFunctionBegin; 2713 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2714 if (u) PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2715 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 */ 2716 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2717 if (!ts->vec_sol || u == ts->vec_sol) { 2718 ierr = VecDuplicate(u,&solution);CHKERRQ(ierr); 2719 ierr = TSSetSolution(ts,solution);CHKERRQ(ierr); 2720 ierr = VecDestroy(&solution);CHKERRQ(ierr); /* grant ownership */ 2721 } 2722 ierr = VecCopy(u,ts->vec_sol);CHKERRQ(ierr); 2723 } else if (u) { 2724 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 2725 } 2726 ierr = TSSetUp(ts);CHKERRQ(ierr); 2727 /* reset time step and iteration counters */ 2728 ts->steps = 0; 2729 ts->ksp_its = 0; 2730 ts->snes_its = 0; 2731 ts->num_snes_failures = 0; 2732 ts->reject = 0; 2733 ts->reason = TS_CONVERGED_ITERATING; 2734 2735 ierr = TSViewFromOptions(ts,NULL,"-ts_view_pre");CHKERRQ(ierr); 2736 2737 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 2738 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 2739 ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr); 2740 ts->solvetime = ts->ptime; 2741 } else { 2742 /* steps the requested number of timesteps. */ 2743 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 2744 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 2745 while (!ts->reason) { 2746 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2747 ierr = TSStep(ts);CHKERRQ(ierr); 2748 ierr = TSPostStep(ts);CHKERRQ(ierr); 2749 } 2750 if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) { 2751 ierr = TSInterpolate(ts,ts->max_time,u);CHKERRQ(ierr); 2752 ts->solvetime = ts->max_time; 2753 solution = u; 2754 } else { 2755 if (u) {ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr);} 2756 ts->solvetime = ts->ptime; 2757 solution = ts->vec_sol; 2758 } 2759 ierr = TSMonitor(ts,ts->steps,ts->solvetime,solution);CHKERRQ(ierr); 2760 } 2761 ierr = TSViewFromOptions(ts,NULL,"-ts_view");CHKERRQ(ierr); 2762 ierr = PetscObjectSAWsBlock((PetscObject)ts);CHKERRQ(ierr); 2763 PetscFunctionReturn(0); 2764 } 2765 2766 #undef __FUNCT__ 2767 #define __FUNCT__ "TSMonitor" 2768 /*@ 2769 TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet() 2770 2771 Collective on TS 2772 2773 Input Parameters: 2774 + ts - time stepping context obtained from TSCreate() 2775 . step - step number that has just completed 2776 . ptime - model time of the state 2777 - u - state at the current model time 2778 2779 Notes: 2780 TSMonitor() is typically used within the time stepping implementations. 2781 Users might call this function when using the TSStep() interface instead of TSSolve(). 2782 2783 Level: advanced 2784 2785 .keywords: TS, timestep 2786 @*/ 2787 PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u) 2788 { 2789 PetscErrorCode ierr; 2790 PetscInt i,n = ts->numbermonitors; 2791 2792 PetscFunctionBegin; 2793 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2794 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 2795 for (i=0; i<n; i++) { 2796 ierr = (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);CHKERRQ(ierr); 2797 } 2798 PetscFunctionReturn(0); 2799 } 2800 2801 /* ------------------------------------------------------------------------*/ 2802 #undef __FUNCT__ 2803 #define __FUNCT__ "TSMonitorLGCtxCreate" 2804 /*@C 2805 TSMonitorLGCtxCreate - Creates a line graph context for use with 2806 TS to monitor the solution process graphically in various ways 2807 2808 Collective on TS 2809 2810 Input Parameters: 2811 + host - the X display to open, or null for the local machine 2812 . label - the title to put in the title bar 2813 . x, y - the screen coordinates of the upper left coordinate of the window 2814 . m, n - the screen width and height in pixels 2815 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time 2816 2817 Output Parameter: 2818 . ctx - the context 2819 2820 Options Database Key: 2821 + -ts_monitor_lg_timestep - automatically sets line graph monitor 2822 . -ts_monitor_lg_solution - 2823 . -ts_monitor_lg_error - 2824 . -ts_monitor_lg_ksp_iterations - 2825 . -ts_monitor_lg_snes_iterations - 2826 - -lg_indicate_data_points <true,false> - indicate the data points (at each time step) on the plot; default is true 2827 2828 Notes: 2829 Use TSMonitorLGCtxDestroy() to destroy. 2830 2831 Level: intermediate 2832 2833 .keywords: TS, monitor, line graph, residual, seealso 2834 2835 .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 2836 2837 @*/ 2838 PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx) 2839 { 2840 PetscDraw win; 2841 PetscErrorCode ierr; 2842 2843 PetscFunctionBegin; 2844 ierr = PetscNew(ctx);CHKERRQ(ierr); 2845 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 2846 ierr = PetscDrawSetFromOptions(win);CHKERRQ(ierr); 2847 ierr = PetscDrawLGCreate(win,1,&(*ctx)->lg);CHKERRQ(ierr); 2848 ierr = PetscLogObjectParent((PetscObject)(*ctx)->lg,(PetscObject)win);CHKERRQ(ierr); 2849 ierr = PetscDrawLGIndicateDataPoints((*ctx)->lg,PETSC_TRUE);CHKERRQ(ierr); 2850 ierr = PetscDrawLGSetFromOptions((*ctx)->lg);CHKERRQ(ierr); 2851 (*ctx)->howoften = howoften; 2852 PetscFunctionReturn(0); 2853 } 2854 2855 #undef __FUNCT__ 2856 #define __FUNCT__ "TSMonitorLGTimeStep" 2857 PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx) 2858 { 2859 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 2860 PetscReal x = ptime,y; 2861 PetscErrorCode ierr; 2862 2863 PetscFunctionBegin; 2864 if (!step) { 2865 PetscDrawAxis axis; 2866 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 2867 ierr = PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time","Time step");CHKERRQ(ierr); 2868 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 2869 ierr = PetscDrawLGIndicateDataPoints(ctx->lg,PETSC_TRUE);CHKERRQ(ierr); 2870 } 2871 ierr = TSGetTimeStep(ts,&y);CHKERRQ(ierr); 2872 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 2873 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 2874 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 2875 } 2876 PetscFunctionReturn(0); 2877 } 2878 2879 #undef __FUNCT__ 2880 #define __FUNCT__ "TSMonitorLGCtxDestroy" 2881 /*@C 2882 TSMonitorLGCtxDestroy - Destroys a line graph context that was created 2883 with TSMonitorLGCtxCreate(). 2884 2885 Collective on TSMonitorLGCtx 2886 2887 Input Parameter: 2888 . ctx - the monitor context 2889 2890 Level: intermediate 2891 2892 .keywords: TS, monitor, line graph, destroy 2893 2894 .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep(); 2895 @*/ 2896 PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx) 2897 { 2898 PetscDraw draw; 2899 PetscErrorCode ierr; 2900 2901 PetscFunctionBegin; 2902 ierr = PetscDrawLGGetDraw((*ctx)->lg,&draw);CHKERRQ(ierr); 2903 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 2904 ierr = PetscDrawLGDestroy(&(*ctx)->lg);CHKERRQ(ierr); 2905 ierr = PetscFree(*ctx);CHKERRQ(ierr); 2906 PetscFunctionReturn(0); 2907 } 2908 2909 #undef __FUNCT__ 2910 #define __FUNCT__ "TSGetTime" 2911 /*@ 2912 TSGetTime - Gets the time of the most recently completed step. 2913 2914 Not Collective 2915 2916 Input Parameter: 2917 . ts - the TS context obtained from TSCreate() 2918 2919 Output Parameter: 2920 . t - the current time 2921 2922 Level: beginner 2923 2924 Note: 2925 When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(), 2926 TSSetPreStage(), TSSetPostStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated. 2927 2928 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 2929 2930 .keywords: TS, get, time 2931 @*/ 2932 PetscErrorCode TSGetTime(TS ts,PetscReal *t) 2933 { 2934 PetscFunctionBegin; 2935 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2936 PetscValidRealPointer(t,2); 2937 *t = ts->ptime; 2938 PetscFunctionReturn(0); 2939 } 2940 2941 #undef __FUNCT__ 2942 #define __FUNCT__ "TSSetTime" 2943 /*@ 2944 TSSetTime - Allows one to reset the time. 2945 2946 Logically Collective on TS 2947 2948 Input Parameters: 2949 + ts - the TS context obtained from TSCreate() 2950 - time - the time 2951 2952 Level: intermediate 2953 2954 .seealso: TSGetTime(), TSSetDuration() 2955 2956 .keywords: TS, set, time 2957 @*/ 2958 PetscErrorCode TSSetTime(TS ts, PetscReal t) 2959 { 2960 PetscFunctionBegin; 2961 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2962 PetscValidLogicalCollectiveReal(ts,t,2); 2963 ts->ptime = t; 2964 PetscFunctionReturn(0); 2965 } 2966 2967 #undef __FUNCT__ 2968 #define __FUNCT__ "TSSetOptionsPrefix" 2969 /*@C 2970 TSSetOptionsPrefix - Sets the prefix used for searching for all 2971 TS options in the database. 2972 2973 Logically Collective on TS 2974 2975 Input Parameter: 2976 + ts - The TS context 2977 - prefix - The prefix to prepend to all option names 2978 2979 Notes: 2980 A hyphen (-) must NOT be given at the beginning of the prefix name. 2981 The first character of all runtime options is AUTOMATICALLY the 2982 hyphen. 2983 2984 Level: advanced 2985 2986 .keywords: TS, set, options, prefix, database 2987 2988 .seealso: TSSetFromOptions() 2989 2990 @*/ 2991 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 2992 { 2993 PetscErrorCode ierr; 2994 SNES snes; 2995 2996 PetscFunctionBegin; 2997 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2998 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2999 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3000 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 3001 PetscFunctionReturn(0); 3002 } 3003 3004 3005 #undef __FUNCT__ 3006 #define __FUNCT__ "TSAppendOptionsPrefix" 3007 /*@C 3008 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 3009 TS options in the database. 3010 3011 Logically Collective on TS 3012 3013 Input Parameter: 3014 + ts - The TS context 3015 - prefix - The prefix to prepend to all option names 3016 3017 Notes: 3018 A hyphen (-) must NOT be given at the beginning of the prefix name. 3019 The first character of all runtime options is AUTOMATICALLY the 3020 hyphen. 3021 3022 Level: advanced 3023 3024 .keywords: TS, append, options, prefix, database 3025 3026 .seealso: TSGetOptionsPrefix() 3027 3028 @*/ 3029 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 3030 { 3031 PetscErrorCode ierr; 3032 SNES snes; 3033 3034 PetscFunctionBegin; 3035 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3036 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3037 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3038 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 3039 PetscFunctionReturn(0); 3040 } 3041 3042 #undef __FUNCT__ 3043 #define __FUNCT__ "TSGetOptionsPrefix" 3044 /*@C 3045 TSGetOptionsPrefix - Sets the prefix used for searching for all 3046 TS options in the database. 3047 3048 Not Collective 3049 3050 Input Parameter: 3051 . ts - The TS context 3052 3053 Output Parameter: 3054 . prefix - A pointer to the prefix string used 3055 3056 Notes: On the fortran side, the user should pass in a string 'prifix' of 3057 sufficient length to hold the prefix. 3058 3059 Level: intermediate 3060 3061 .keywords: TS, get, options, prefix, database 3062 3063 .seealso: TSAppendOptionsPrefix() 3064 @*/ 3065 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 3066 { 3067 PetscErrorCode ierr; 3068 3069 PetscFunctionBegin; 3070 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3071 PetscValidPointer(prefix,2); 3072 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3073 PetscFunctionReturn(0); 3074 } 3075 3076 #undef __FUNCT__ 3077 #define __FUNCT__ "TSGetRHSJacobian" 3078 /*@C 3079 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 3080 3081 Not Collective, but parallel objects are returned if TS is parallel 3082 3083 Input Parameter: 3084 . ts - The TS context obtained from TSCreate() 3085 3086 Output Parameters: 3087 + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or NULL) 3088 . Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat (or NULL) 3089 . func - Function to compute the Jacobian of the RHS (or NULL) 3090 - ctx - User-defined context for Jacobian evaluation routine (or NULL) 3091 3092 Notes: You can pass in NULL for any return argument you do not need. 3093 3094 Level: intermediate 3095 3096 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3097 3098 .keywords: TS, timestep, get, matrix, Jacobian 3099 @*/ 3100 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx) 3101 { 3102 PetscErrorCode ierr; 3103 SNES snes; 3104 DM dm; 3105 3106 PetscFunctionBegin; 3107 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3108 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3109 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3110 ierr = DMTSGetRHSJacobian(dm,func,ctx);CHKERRQ(ierr); 3111 PetscFunctionReturn(0); 3112 } 3113 3114 #undef __FUNCT__ 3115 #define __FUNCT__ "TSGetIJacobian" 3116 /*@C 3117 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 3118 3119 Not Collective, but parallel objects are returned if TS is parallel 3120 3121 Input Parameter: 3122 . ts - The TS context obtained from TSCreate() 3123 3124 Output Parameters: 3125 + Amat - The (approximate) Jacobian of F(t,U,U_t) 3126 . Pmat - The matrix from which the preconditioner is constructed, often the same as Amat 3127 . f - The function to compute the matrices 3128 - ctx - User-defined context for Jacobian evaluation routine 3129 3130 Notes: You can pass in NULL for any return argument you do not need. 3131 3132 Level: advanced 3133 3134 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3135 3136 .keywords: TS, timestep, get, matrix, Jacobian 3137 @*/ 3138 PetscErrorCode TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx) 3139 { 3140 PetscErrorCode ierr; 3141 SNES snes; 3142 DM dm; 3143 3144 PetscFunctionBegin; 3145 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3146 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 3147 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3148 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3149 ierr = DMTSGetIJacobian(dm,f,ctx);CHKERRQ(ierr); 3150 PetscFunctionReturn(0); 3151 } 3152 3153 3154 #undef __FUNCT__ 3155 #define __FUNCT__ "TSMonitorDrawSolution" 3156 /*@C 3157 TSMonitorDrawSolution - Monitors progress of the TS solvers by calling 3158 VecView() for the solution at each timestep 3159 3160 Collective on TS 3161 3162 Input Parameters: 3163 + ts - the TS context 3164 . step - current time-step 3165 . ptime - current time 3166 - dummy - either a viewer or NULL 3167 3168 Options Database: 3169 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3170 3171 Notes: the initial solution and current solution are not displayed with a common axis scaling so generally the option -ts_monitor_draw_solution_initial 3172 will look bad 3173 3174 Level: intermediate 3175 3176 .keywords: TS, vector, monitor, view 3177 3178 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3179 @*/ 3180 PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3181 { 3182 PetscErrorCode ierr; 3183 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3184 PetscDraw draw; 3185 3186 PetscFunctionBegin; 3187 if (!step && ictx->showinitial) { 3188 if (!ictx->initialsolution) { 3189 ierr = VecDuplicate(u,&ictx->initialsolution);CHKERRQ(ierr); 3190 } 3191 ierr = VecCopy(u,ictx->initialsolution);CHKERRQ(ierr); 3192 } 3193 if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 3194 3195 if (ictx->showinitial) { 3196 PetscReal pause; 3197 ierr = PetscViewerDrawGetPause(ictx->viewer,&pause);CHKERRQ(ierr); 3198 ierr = PetscViewerDrawSetPause(ictx->viewer,0.0);CHKERRQ(ierr); 3199 ierr = VecView(ictx->initialsolution,ictx->viewer);CHKERRQ(ierr); 3200 ierr = PetscViewerDrawSetPause(ictx->viewer,pause);CHKERRQ(ierr); 3201 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);CHKERRQ(ierr); 3202 } 3203 ierr = VecView(u,ictx->viewer);CHKERRQ(ierr); 3204 if (ictx->showtimestepandtime) { 3205 PetscReal xl,yl,xr,yr,tw,w,h; 3206 char time[32]; 3207 size_t len; 3208 3209 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3210 ierr = PetscSNPrintf(time,32,"Timestep %d Time %f",(int)step,(double)ptime);CHKERRQ(ierr); 3211 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3212 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3213 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3214 w = xl + .5*(xr - xl) - .5*len*tw; 3215 h = yl + .95*(yr - yl); 3216 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3217 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3218 } 3219 3220 if (ictx->showinitial) { 3221 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);CHKERRQ(ierr); 3222 } 3223 PetscFunctionReturn(0); 3224 } 3225 3226 #undef __FUNCT__ 3227 #define __FUNCT__ "TSMonitorDrawSolutionPhase" 3228 /*@C 3229 TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram 3230 3231 Collective on TS 3232 3233 Input Parameters: 3234 + ts - the TS context 3235 . step - current time-step 3236 . ptime - current time 3237 - dummy - either a viewer or NULL 3238 3239 Level: intermediate 3240 3241 .keywords: TS, vector, monitor, view 3242 3243 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3244 @*/ 3245 PetscErrorCode TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3246 { 3247 PetscErrorCode ierr; 3248 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3249 PetscDraw draw; 3250 MPI_Comm comm; 3251 PetscInt n; 3252 PetscMPIInt size; 3253 PetscReal xl,yl,xr,yr,tw,w,h; 3254 char time[32]; 3255 size_t len; 3256 const PetscScalar *U; 3257 3258 PetscFunctionBegin; 3259 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 3260 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3261 if (size != 1) SETERRQ(comm,PETSC_ERR_SUP,"Only allowed for sequential runs"); 3262 ierr = VecGetSize(u,&n);CHKERRQ(ierr); 3263 if (n != 2) SETERRQ(comm,PETSC_ERR_SUP,"Only for ODEs with two unknowns"); 3264 3265 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3266 3267 ierr = VecGetArrayRead(u,&U);CHKERRQ(ierr); 3268 ierr = PetscDrawAxisGetLimits(ictx->axis,&xl,&xr,&yl,&yr);CHKERRQ(ierr); 3269 if ((PetscRealPart(U[0]) < xl) || (PetscRealPart(U[1]) < yl) || (PetscRealPart(U[0]) > xr) || (PetscRealPart(U[1]) > yr)) { 3270 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3271 PetscFunctionReturn(0); 3272 } 3273 if (!step) ictx->color++; 3274 ierr = PetscDrawPoint(draw,PetscRealPart(U[0]),PetscRealPart(U[1]),ictx->color);CHKERRQ(ierr); 3275 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3276 3277 if (ictx->showtimestepandtime) { 3278 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3279 ierr = PetscSNPrintf(time,32,"Timestep %d Time %f",(int)step,(double)ptime);CHKERRQ(ierr); 3280 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3281 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3282 w = xl + .5*(xr - xl) - .5*len*tw; 3283 h = yl + .95*(yr - yl); 3284 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3285 } 3286 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3287 PetscFunctionReturn(0); 3288 } 3289 3290 3291 #undef __FUNCT__ 3292 #define __FUNCT__ "TSMonitorDrawCtxDestroy" 3293 /*@C 3294 TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution() 3295 3296 Collective on TS 3297 3298 Input Parameters: 3299 . ctx - the monitor context 3300 3301 Level: intermediate 3302 3303 .keywords: TS, vector, monitor, view 3304 3305 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError() 3306 @*/ 3307 PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx) 3308 { 3309 PetscErrorCode ierr; 3310 3311 PetscFunctionBegin; 3312 ierr = PetscDrawAxisDestroy(&(*ictx)->axis);CHKERRQ(ierr); 3313 ierr = PetscViewerDestroy(&(*ictx)->viewer);CHKERRQ(ierr); 3314 ierr = VecDestroy(&(*ictx)->initialsolution);CHKERRQ(ierr); 3315 ierr = PetscFree(*ictx);CHKERRQ(ierr); 3316 PetscFunctionReturn(0); 3317 } 3318 3319 #undef __FUNCT__ 3320 #define __FUNCT__ "TSMonitorDrawCtxCreate" 3321 /*@C 3322 TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx 3323 3324 Collective on TS 3325 3326 Input Parameter: 3327 . ts - time-step context 3328 3329 Output Patameter: 3330 . ctx - the monitor context 3331 3332 Options Database: 3333 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3334 3335 Level: intermediate 3336 3337 .keywords: TS, vector, monitor, view 3338 3339 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx() 3340 @*/ 3341 PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx) 3342 { 3343 PetscErrorCode ierr; 3344 3345 PetscFunctionBegin; 3346 ierr = PetscNew(ctx);CHKERRQ(ierr); 3347 ierr = PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);CHKERRQ(ierr); 3348 ierr = PetscViewerSetFromOptions((*ctx)->viewer);CHKERRQ(ierr); 3349 3350 (*ctx)->howoften = howoften; 3351 (*ctx)->showinitial = PETSC_FALSE; 3352 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);CHKERRQ(ierr); 3353 3354 (*ctx)->showtimestepandtime = PETSC_FALSE; 3355 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);CHKERRQ(ierr); 3356 (*ctx)->color = PETSC_DRAW_WHITE; 3357 PetscFunctionReturn(0); 3358 } 3359 3360 #undef __FUNCT__ 3361 #define __FUNCT__ "TSMonitorDrawError" 3362 /*@C 3363 TSMonitorDrawError - Monitors progress of the TS solvers by calling 3364 VecView() for the error at each timestep 3365 3366 Collective on TS 3367 3368 Input Parameters: 3369 + ts - the TS context 3370 . step - current time-step 3371 . ptime - current time 3372 - dummy - either a viewer or NULL 3373 3374 Level: intermediate 3375 3376 .keywords: TS, vector, monitor, view 3377 3378 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3379 @*/ 3380 PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3381 { 3382 PetscErrorCode ierr; 3383 TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy; 3384 PetscViewer viewer = ctx->viewer; 3385 Vec work; 3386 3387 PetscFunctionBegin; 3388 if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 3389 ierr = VecDuplicate(u,&work);CHKERRQ(ierr); 3390 ierr = TSComputeSolutionFunction(ts,ptime,work);CHKERRQ(ierr); 3391 ierr = VecAXPY(work,-1.0,u);CHKERRQ(ierr); 3392 ierr = VecView(work,viewer);CHKERRQ(ierr); 3393 ierr = VecDestroy(&work);CHKERRQ(ierr); 3394 PetscFunctionReturn(0); 3395 } 3396 3397 #include <petsc-private/dmimpl.h> 3398 #undef __FUNCT__ 3399 #define __FUNCT__ "TSSetDM" 3400 /*@ 3401 TSSetDM - Sets the DM that may be used by some preconditioners 3402 3403 Logically Collective on TS and DM 3404 3405 Input Parameters: 3406 + ts - the preconditioner context 3407 - dm - the dm 3408 3409 Level: intermediate 3410 3411 3412 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 3413 @*/ 3414 PetscErrorCode TSSetDM(TS ts,DM dm) 3415 { 3416 PetscErrorCode ierr; 3417 SNES snes; 3418 DMTS tsdm; 3419 3420 PetscFunctionBegin; 3421 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3422 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 3423 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 3424 if (ts->dm->dmts && !dm->dmts) { 3425 ierr = DMCopyDMTS(ts->dm,dm);CHKERRQ(ierr); 3426 ierr = DMGetDMTS(ts->dm,&tsdm);CHKERRQ(ierr); 3427 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 3428 tsdm->originaldm = dm; 3429 } 3430 } 3431 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 3432 } 3433 ts->dm = dm; 3434 3435 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3436 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 3437 PetscFunctionReturn(0); 3438 } 3439 3440 #undef __FUNCT__ 3441 #define __FUNCT__ "TSGetDM" 3442 /*@ 3443 TSGetDM - Gets the DM that may be used by some preconditioners 3444 3445 Not Collective 3446 3447 Input Parameter: 3448 . ts - the preconditioner context 3449 3450 Output Parameter: 3451 . dm - the dm 3452 3453 Level: intermediate 3454 3455 3456 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 3457 @*/ 3458 PetscErrorCode TSGetDM(TS ts,DM *dm) 3459 { 3460 PetscErrorCode ierr; 3461 3462 PetscFunctionBegin; 3463 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3464 if (!ts->dm) { 3465 ierr = DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);CHKERRQ(ierr); 3466 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 3467 } 3468 *dm = ts->dm; 3469 PetscFunctionReturn(0); 3470 } 3471 3472 #undef __FUNCT__ 3473 #define __FUNCT__ "SNESTSFormFunction" 3474 /*@ 3475 SNESTSFormFunction - Function to evaluate nonlinear residual 3476 3477 Logically Collective on SNES 3478 3479 Input Parameter: 3480 + snes - nonlinear solver 3481 . U - the current state at which to evaluate the residual 3482 - ctx - user context, must be a TS 3483 3484 Output Parameter: 3485 . F - the nonlinear residual 3486 3487 Notes: 3488 This function is not normally called by users and is automatically registered with the SNES used by TS. 3489 It is most frequently passed to MatFDColoringSetFunction(). 3490 3491 Level: advanced 3492 3493 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 3494 @*/ 3495 PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx) 3496 { 3497 TS ts = (TS)ctx; 3498 PetscErrorCode ierr; 3499 3500 PetscFunctionBegin; 3501 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3502 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3503 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 3504 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 3505 ierr = (ts->ops->snesfunction)(snes,U,F,ts);CHKERRQ(ierr); 3506 PetscFunctionReturn(0); 3507 } 3508 3509 #undef __FUNCT__ 3510 #define __FUNCT__ "SNESTSFormJacobian" 3511 /*@ 3512 SNESTSFormJacobian - Function to evaluate the Jacobian 3513 3514 Collective on SNES 3515 3516 Input Parameter: 3517 + snes - nonlinear solver 3518 . U - the current state at which to evaluate the residual 3519 - ctx - user context, must be a TS 3520 3521 Output Parameter: 3522 + A - the Jacobian 3523 . B - the preconditioning matrix (may be the same as A) 3524 - flag - indicates any structure change in the matrix 3525 3526 Notes: 3527 This function is not normally called by users and is automatically registered with the SNES used by TS. 3528 3529 Level: developer 3530 3531 .seealso: SNESSetJacobian() 3532 @*/ 3533 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat A,Mat B,void *ctx) 3534 { 3535 TS ts = (TS)ctx; 3536 PetscErrorCode ierr; 3537 3538 PetscFunctionBegin; 3539 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3540 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3541 PetscValidPointer(A,3); 3542 PetscValidHeaderSpecific(A,MAT_CLASSID,3); 3543 PetscValidPointer(B,4); 3544 PetscValidHeaderSpecific(B,MAT_CLASSID,4); 3545 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 3546 ierr = (ts->ops->snesjacobian)(snes,U,A,B,ts);CHKERRQ(ierr); 3547 PetscFunctionReturn(0); 3548 } 3549 3550 #undef __FUNCT__ 3551 #define __FUNCT__ "TSComputeRHSFunctionLinear" 3552 /*@C 3553 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 3554 3555 Collective on TS 3556 3557 Input Arguments: 3558 + ts - time stepping context 3559 . t - time at which to evaluate 3560 . U - state at which to evaluate 3561 - ctx - context 3562 3563 Output Arguments: 3564 . F - right hand side 3565 3566 Level: intermediate 3567 3568 Notes: 3569 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 3570 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 3571 3572 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 3573 @*/ 3574 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 3575 { 3576 PetscErrorCode ierr; 3577 Mat Arhs,Brhs; 3578 3579 PetscFunctionBegin; 3580 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 3581 ierr = TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);CHKERRQ(ierr); 3582 ierr = MatMult(Arhs,U,F);CHKERRQ(ierr); 3583 PetscFunctionReturn(0); 3584 } 3585 3586 #undef __FUNCT__ 3587 #define __FUNCT__ "TSComputeRHSJacobianConstant" 3588 /*@C 3589 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 3590 3591 Collective on TS 3592 3593 Input Arguments: 3594 + ts - time stepping context 3595 . t - time at which to evaluate 3596 . U - state at which to evaluate 3597 - ctx - context 3598 3599 Output Arguments: 3600 + A - pointer to operator 3601 . B - pointer to preconditioning matrix 3602 - flg - matrix structure flag 3603 3604 Level: intermediate 3605 3606 Notes: 3607 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 3608 3609 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 3610 @*/ 3611 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx) 3612 { 3613 PetscFunctionBegin; 3614 PetscFunctionReturn(0); 3615 } 3616 3617 #undef __FUNCT__ 3618 #define __FUNCT__ "TSComputeIFunctionLinear" 3619 /*@C 3620 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 3621 3622 Collective on TS 3623 3624 Input Arguments: 3625 + ts - time stepping context 3626 . t - time at which to evaluate 3627 . U - state at which to evaluate 3628 . Udot - time derivative of state vector 3629 - ctx - context 3630 3631 Output Arguments: 3632 . F - left hand side 3633 3634 Level: intermediate 3635 3636 Notes: 3637 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 3638 user is required to write their own TSComputeIFunction. 3639 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 3640 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 3641 3642 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 3643 @*/ 3644 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 3645 { 3646 PetscErrorCode ierr; 3647 Mat A,B; 3648 3649 PetscFunctionBegin; 3650 ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr); 3651 ierr = TSComputeIJacobian(ts,t,U,Udot,1.0,A,B,PETSC_TRUE);CHKERRQ(ierr); 3652 ierr = MatMult(A,Udot,F);CHKERRQ(ierr); 3653 PetscFunctionReturn(0); 3654 } 3655 3656 #undef __FUNCT__ 3657 #define __FUNCT__ "TSComputeIJacobianConstant" 3658 /*@C 3659 TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE 3660 3661 Collective on TS 3662 3663 Input Arguments: 3664 + ts - time stepping context 3665 . t - time at which to evaluate 3666 . U - state at which to evaluate 3667 . Udot - time derivative of state vector 3668 . shift - shift to apply 3669 - ctx - context 3670 3671 Output Arguments: 3672 + A - pointer to operator 3673 . B - pointer to preconditioning matrix 3674 - flg - matrix structure flag 3675 3676 Level: advanced 3677 3678 Notes: 3679 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 3680 3681 It is only appropriate for problems of the form 3682 3683 $ M Udot = F(U,t) 3684 3685 where M is constant and F is non-stiff. The user must pass M to TSSetIJacobian(). The current implementation only 3686 works with IMEX time integration methods such as TSROSW and TSARKIMEX, since there is no support for de-constructing 3687 an implicit operator of the form 3688 3689 $ shift*M + J 3690 3691 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 3692 a copy of M or reassemble it when requested. 3693 3694 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 3695 @*/ 3696 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,void *ctx) 3697 { 3698 PetscErrorCode ierr; 3699 3700 PetscFunctionBegin; 3701 ierr = MatScale(A, shift / ts->ijacobian.shift);CHKERRQ(ierr); 3702 ts->ijacobian.shift = shift; 3703 PetscFunctionReturn(0); 3704 } 3705 3706 #undef __FUNCT__ 3707 #define __FUNCT__ "TSGetEquationType" 3708 /*@ 3709 TSGetEquationType - Gets the type of the equation that TS is solving. 3710 3711 Not Collective 3712 3713 Input Parameter: 3714 . ts - the TS context 3715 3716 Output Parameter: 3717 . equation_type - see TSEquationType 3718 3719 Level: beginner 3720 3721 .keywords: TS, equation type 3722 3723 .seealso: TSSetEquationType(), TSEquationType 3724 @*/ 3725 PetscErrorCode TSGetEquationType(TS ts,TSEquationType *equation_type) 3726 { 3727 PetscFunctionBegin; 3728 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3729 PetscValidPointer(equation_type,2); 3730 *equation_type = ts->equation_type; 3731 PetscFunctionReturn(0); 3732 } 3733 3734 #undef __FUNCT__ 3735 #define __FUNCT__ "TSSetEquationType" 3736 /*@ 3737 TSSetEquationType - Sets the type of the equation that TS is solving. 3738 3739 Not Collective 3740 3741 Input Parameter: 3742 + ts - the TS context 3743 . equation_type - see TSEquationType 3744 3745 Level: advanced 3746 3747 .keywords: TS, equation type 3748 3749 .seealso: TSGetEquationType(), TSEquationType 3750 @*/ 3751 PetscErrorCode TSSetEquationType(TS ts,TSEquationType equation_type) 3752 { 3753 PetscFunctionBegin; 3754 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3755 ts->equation_type = equation_type; 3756 PetscFunctionReturn(0); 3757 } 3758 3759 #undef __FUNCT__ 3760 #define __FUNCT__ "TSGetConvergedReason" 3761 /*@ 3762 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 3763 3764 Not Collective 3765 3766 Input Parameter: 3767 . ts - the TS context 3768 3769 Output Parameter: 3770 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3771 manual pages for the individual convergence tests for complete lists 3772 3773 Level: beginner 3774 3775 Notes: 3776 Can only be called after the call to TSSolve() is complete. 3777 3778 .keywords: TS, nonlinear, set, convergence, test 3779 3780 .seealso: TSSetConvergenceTest(), TSConvergedReason 3781 @*/ 3782 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 3783 { 3784 PetscFunctionBegin; 3785 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3786 PetscValidPointer(reason,2); 3787 *reason = ts->reason; 3788 PetscFunctionReturn(0); 3789 } 3790 3791 #undef __FUNCT__ 3792 #define __FUNCT__ "TSSetConvergedReason" 3793 /*@ 3794 TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve. 3795 3796 Not Collective 3797 3798 Input Parameter: 3799 + ts - the TS context 3800 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3801 manual pages for the individual convergence tests for complete lists 3802 3803 Level: advanced 3804 3805 Notes: 3806 Can only be called during TSSolve() is active. 3807 3808 .keywords: TS, nonlinear, set, convergence, test 3809 3810 .seealso: TSConvergedReason 3811 @*/ 3812 PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason) 3813 { 3814 PetscFunctionBegin; 3815 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3816 ts->reason = reason; 3817 PetscFunctionReturn(0); 3818 } 3819 3820 #undef __FUNCT__ 3821 #define __FUNCT__ "TSGetSolveTime" 3822 /*@ 3823 TSGetSolveTime - Gets the time after a call to TSSolve() 3824 3825 Not Collective 3826 3827 Input Parameter: 3828 . ts - the TS context 3829 3830 Output Parameter: 3831 . ftime - the final time. This time should correspond to the final time set with TSSetDuration() 3832 3833 Level: beginner 3834 3835 Notes: 3836 Can only be called after the call to TSSolve() is complete. 3837 3838 .keywords: TS, nonlinear, set, convergence, test 3839 3840 .seealso: TSSetConvergenceTest(), TSConvergedReason 3841 @*/ 3842 PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime) 3843 { 3844 PetscFunctionBegin; 3845 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3846 PetscValidPointer(ftime,2); 3847 *ftime = ts->solvetime; 3848 PetscFunctionReturn(0); 3849 } 3850 3851 #undef __FUNCT__ 3852 #define __FUNCT__ "TSGetSNESIterations" 3853 /*@ 3854 TSGetSNESIterations - Gets the total number of nonlinear iterations 3855 used by the time integrator. 3856 3857 Not Collective 3858 3859 Input Parameter: 3860 . ts - TS context 3861 3862 Output Parameter: 3863 . nits - number of nonlinear iterations 3864 3865 Notes: 3866 This counter is reset to zero for each successive call to TSSolve(). 3867 3868 Level: intermediate 3869 3870 .keywords: TS, get, number, nonlinear, iterations 3871 3872 .seealso: TSGetKSPIterations() 3873 @*/ 3874 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 3875 { 3876 PetscFunctionBegin; 3877 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3878 PetscValidIntPointer(nits,2); 3879 *nits = ts->snes_its; 3880 PetscFunctionReturn(0); 3881 } 3882 3883 #undef __FUNCT__ 3884 #define __FUNCT__ "TSGetKSPIterations" 3885 /*@ 3886 TSGetKSPIterations - Gets the total number of linear iterations 3887 used by the time integrator. 3888 3889 Not Collective 3890 3891 Input Parameter: 3892 . ts - TS context 3893 3894 Output Parameter: 3895 . lits - number of linear iterations 3896 3897 Notes: 3898 This counter is reset to zero for each successive call to TSSolve(). 3899 3900 Level: intermediate 3901 3902 .keywords: TS, get, number, linear, iterations 3903 3904 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 3905 @*/ 3906 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 3907 { 3908 PetscFunctionBegin; 3909 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3910 PetscValidIntPointer(lits,2); 3911 *lits = ts->ksp_its; 3912 PetscFunctionReturn(0); 3913 } 3914 3915 #undef __FUNCT__ 3916 #define __FUNCT__ "TSGetStepRejections" 3917 /*@ 3918 TSGetStepRejections - Gets the total number of rejected steps. 3919 3920 Not Collective 3921 3922 Input Parameter: 3923 . ts - TS context 3924 3925 Output Parameter: 3926 . rejects - number of steps rejected 3927 3928 Notes: 3929 This counter is reset to zero for each successive call to TSSolve(). 3930 3931 Level: intermediate 3932 3933 .keywords: TS, get, number 3934 3935 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 3936 @*/ 3937 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 3938 { 3939 PetscFunctionBegin; 3940 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3941 PetscValidIntPointer(rejects,2); 3942 *rejects = ts->reject; 3943 PetscFunctionReturn(0); 3944 } 3945 3946 #undef __FUNCT__ 3947 #define __FUNCT__ "TSGetSNESFailures" 3948 /*@ 3949 TSGetSNESFailures - Gets the total number of failed SNES solves 3950 3951 Not Collective 3952 3953 Input Parameter: 3954 . ts - TS context 3955 3956 Output Parameter: 3957 . fails - number of failed nonlinear solves 3958 3959 Notes: 3960 This counter is reset to zero for each successive call to TSSolve(). 3961 3962 Level: intermediate 3963 3964 .keywords: TS, get, number 3965 3966 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 3967 @*/ 3968 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 3969 { 3970 PetscFunctionBegin; 3971 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3972 PetscValidIntPointer(fails,2); 3973 *fails = ts->num_snes_failures; 3974 PetscFunctionReturn(0); 3975 } 3976 3977 #undef __FUNCT__ 3978 #define __FUNCT__ "TSSetMaxStepRejections" 3979 /*@ 3980 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 3981 3982 Not Collective 3983 3984 Input Parameter: 3985 + ts - TS context 3986 - rejects - maximum number of rejected steps, pass -1 for unlimited 3987 3988 Notes: 3989 The counter is reset to zero for each step 3990 3991 Options Database Key: 3992 . -ts_max_reject - Maximum number of step rejections before a step fails 3993 3994 Level: intermediate 3995 3996 .keywords: TS, set, maximum, number 3997 3998 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 3999 @*/ 4000 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 4001 { 4002 PetscFunctionBegin; 4003 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4004 ts->max_reject = rejects; 4005 PetscFunctionReturn(0); 4006 } 4007 4008 #undef __FUNCT__ 4009 #define __FUNCT__ "TSSetMaxSNESFailures" 4010 /*@ 4011 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 4012 4013 Not Collective 4014 4015 Input Parameter: 4016 + ts - TS context 4017 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 4018 4019 Notes: 4020 The counter is reset to zero for each successive call to TSSolve(). 4021 4022 Options Database Key: 4023 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 4024 4025 Level: intermediate 4026 4027 .keywords: TS, set, maximum, number 4028 4029 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 4030 @*/ 4031 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 4032 { 4033 PetscFunctionBegin; 4034 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4035 ts->max_snes_failures = fails; 4036 PetscFunctionReturn(0); 4037 } 4038 4039 #undef __FUNCT__ 4040 #define __FUNCT__ "TSSetErrorIfStepFails" 4041 /*@ 4042 TSSetErrorIfStepFails - Error if no step succeeds 4043 4044 Not Collective 4045 4046 Input Parameter: 4047 + ts - TS context 4048 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 4049 4050 Options Database Key: 4051 . -ts_error_if_step_fails - Error if no step succeeds 4052 4053 Level: intermediate 4054 4055 .keywords: TS, set, error 4056 4057 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 4058 @*/ 4059 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 4060 { 4061 PetscFunctionBegin; 4062 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4063 ts->errorifstepfailed = err; 4064 PetscFunctionReturn(0); 4065 } 4066 4067 #undef __FUNCT__ 4068 #define __FUNCT__ "TSMonitorSolutionBinary" 4069 /*@C 4070 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 4071 4072 Collective on TS 4073 4074 Input Parameters: 4075 + ts - the TS context 4076 . step - current time-step 4077 . ptime - current time 4078 . u - current state 4079 - viewer - binary viewer 4080 4081 Level: intermediate 4082 4083 .keywords: TS, vector, monitor, view 4084 4085 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4086 @*/ 4087 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec u,void *viewer) 4088 { 4089 PetscErrorCode ierr; 4090 PetscViewer v = (PetscViewer)viewer; 4091 4092 PetscFunctionBegin; 4093 ierr = VecView(u,v);CHKERRQ(ierr); 4094 PetscFunctionReturn(0); 4095 } 4096 4097 #undef __FUNCT__ 4098 #define __FUNCT__ "TSMonitorSolutionVTK" 4099 /*@C 4100 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 4101 4102 Collective on TS 4103 4104 Input Parameters: 4105 + ts - the TS context 4106 . step - current time-step 4107 . ptime - current time 4108 . u - current state 4109 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4110 4111 Level: intermediate 4112 4113 Notes: 4114 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. 4115 These are named according to the file name template. 4116 4117 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 4118 4119 .keywords: TS, vector, monitor, view 4120 4121 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4122 @*/ 4123 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate) 4124 { 4125 PetscErrorCode ierr; 4126 char filename[PETSC_MAX_PATH_LEN]; 4127 PetscViewer viewer; 4128 4129 PetscFunctionBegin; 4130 ierr = PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);CHKERRQ(ierr); 4131 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 4132 ierr = VecView(u,viewer);CHKERRQ(ierr); 4133 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4134 PetscFunctionReturn(0); 4135 } 4136 4137 #undef __FUNCT__ 4138 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 4139 /*@C 4140 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 4141 4142 Collective on TS 4143 4144 Input Parameters: 4145 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4146 4147 Level: intermediate 4148 4149 Note: 4150 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 4151 4152 .keywords: TS, vector, monitor, view 4153 4154 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 4155 @*/ 4156 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 4157 { 4158 PetscErrorCode ierr; 4159 4160 PetscFunctionBegin; 4161 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 4162 PetscFunctionReturn(0); 4163 } 4164 4165 #undef __FUNCT__ 4166 #define __FUNCT__ "TSGetAdapt" 4167 /*@ 4168 TSGetAdapt - Get the adaptive controller context for the current method 4169 4170 Collective on TS if controller has not been created yet 4171 4172 Input Arguments: 4173 . ts - time stepping context 4174 4175 Output Arguments: 4176 . adapt - adaptive controller 4177 4178 Level: intermediate 4179 4180 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 4181 @*/ 4182 PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt) 4183 { 4184 PetscErrorCode ierr; 4185 4186 PetscFunctionBegin; 4187 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4188 PetscValidPointer(adapt,2); 4189 if (!ts->adapt) { 4190 ierr = TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);CHKERRQ(ierr); 4191 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->adapt);CHKERRQ(ierr); 4192 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 4193 } 4194 *adapt = ts->adapt; 4195 PetscFunctionReturn(0); 4196 } 4197 4198 #undef __FUNCT__ 4199 #define __FUNCT__ "TSSetTolerances" 4200 /*@ 4201 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 4202 4203 Logically Collective 4204 4205 Input Arguments: 4206 + ts - time integration context 4207 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 4208 . vatol - vector of absolute tolerances or NULL, used in preference to atol if present 4209 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 4210 - vrtol - vector of relative tolerances or NULL, used in preference to atol if present 4211 4212 Level: beginner 4213 4214 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 4215 @*/ 4216 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 4217 { 4218 PetscErrorCode ierr; 4219 4220 PetscFunctionBegin; 4221 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 4222 if (vatol) { 4223 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 4224 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 4225 4226 ts->vatol = vatol; 4227 } 4228 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 4229 if (vrtol) { 4230 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 4231 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 4232 4233 ts->vrtol = vrtol; 4234 } 4235 PetscFunctionReturn(0); 4236 } 4237 4238 #undef __FUNCT__ 4239 #define __FUNCT__ "TSGetTolerances" 4240 /*@ 4241 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 4242 4243 Logically Collective 4244 4245 Input Arguments: 4246 . ts - time integration context 4247 4248 Output Arguments: 4249 + atol - scalar absolute tolerances, NULL to ignore 4250 . vatol - vector of absolute tolerances, NULL to ignore 4251 . rtol - scalar relative tolerances, NULL to ignore 4252 - vrtol - vector of relative tolerances, NULL to ignore 4253 4254 Level: beginner 4255 4256 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 4257 @*/ 4258 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 4259 { 4260 PetscFunctionBegin; 4261 if (atol) *atol = ts->atol; 4262 if (vatol) *vatol = ts->vatol; 4263 if (rtol) *rtol = ts->rtol; 4264 if (vrtol) *vrtol = ts->vrtol; 4265 PetscFunctionReturn(0); 4266 } 4267 4268 #undef __FUNCT__ 4269 #define __FUNCT__ "TSErrorNormWRMS" 4270 /*@ 4271 TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state 4272 4273 Collective on TS 4274 4275 Input Arguments: 4276 + ts - time stepping context 4277 - Y - state vector to be compared to ts->vec_sol 4278 4279 Output Arguments: 4280 . norm - weighted norm, a value of 1.0 is considered small 4281 4282 Level: developer 4283 4284 .seealso: TSSetTolerances() 4285 @*/ 4286 PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm) 4287 { 4288 PetscErrorCode ierr; 4289 PetscInt i,n,N; 4290 const PetscScalar *u,*y; 4291 Vec U; 4292 PetscReal sum,gsum; 4293 4294 PetscFunctionBegin; 4295 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4296 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 4297 PetscValidPointer(norm,3); 4298 U = ts->vec_sol; 4299 PetscCheckSameTypeAndComm(U,1,Y,2); 4300 if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 4301 4302 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 4303 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 4304 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 4305 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 4306 sum = 0.; 4307 if (ts->vatol && ts->vrtol) { 4308 const PetscScalar *atol,*rtol; 4309 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4310 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4311 for (i=0; i<n; i++) { 4312 PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4313 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4314 } 4315 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4316 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4317 } else if (ts->vatol) { /* vector atol, scalar rtol */ 4318 const PetscScalar *atol; 4319 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4320 for (i=0; i<n; i++) { 4321 PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4322 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4323 } 4324 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4325 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 4326 const PetscScalar *rtol; 4327 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4328 for (i=0; i<n; i++) { 4329 PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4330 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4331 } 4332 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4333 } else { /* scalar atol, scalar rtol */ 4334 for (i=0; i<n; i++) { 4335 PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4336 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4337 } 4338 } 4339 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 4340 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 4341 4342 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4343 *norm = PetscSqrtReal(gsum / N); 4344 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 4345 PetscFunctionReturn(0); 4346 } 4347 4348 #undef __FUNCT__ 4349 #define __FUNCT__ "TSSetCFLTimeLocal" 4350 /*@ 4351 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 4352 4353 Logically Collective on TS 4354 4355 Input Arguments: 4356 + ts - time stepping context 4357 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 4358 4359 Note: 4360 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 4361 4362 Level: intermediate 4363 4364 .seealso: TSGetCFLTime(), TSADAPTCFL 4365 @*/ 4366 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 4367 { 4368 PetscFunctionBegin; 4369 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4370 ts->cfltime_local = cfltime; 4371 ts->cfltime = -1.; 4372 PetscFunctionReturn(0); 4373 } 4374 4375 #undef __FUNCT__ 4376 #define __FUNCT__ "TSGetCFLTime" 4377 /*@ 4378 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 4379 4380 Collective on TS 4381 4382 Input Arguments: 4383 . ts - time stepping context 4384 4385 Output Arguments: 4386 . cfltime - maximum stable time step for forward Euler 4387 4388 Level: advanced 4389 4390 .seealso: TSSetCFLTimeLocal() 4391 @*/ 4392 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 4393 { 4394 PetscErrorCode ierr; 4395 4396 PetscFunctionBegin; 4397 if (ts->cfltime < 0) { 4398 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4399 } 4400 *cfltime = ts->cfltime; 4401 PetscFunctionReturn(0); 4402 } 4403 4404 #undef __FUNCT__ 4405 #define __FUNCT__ "TSVISetVariableBounds" 4406 /*@ 4407 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 4408 4409 Input Parameters: 4410 . ts - the TS context. 4411 . xl - lower bound. 4412 . xu - upper bound. 4413 4414 Notes: 4415 If this routine is not called then the lower and upper bounds are set to 4416 PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp(). 4417 4418 Level: advanced 4419 4420 @*/ 4421 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 4422 { 4423 PetscErrorCode ierr; 4424 SNES snes; 4425 4426 PetscFunctionBegin; 4427 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4428 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 4429 PetscFunctionReturn(0); 4430 } 4431 4432 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4433 #include <mex.h> 4434 4435 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 4436 4437 #undef __FUNCT__ 4438 #define __FUNCT__ "TSComputeFunction_Matlab" 4439 /* 4440 TSComputeFunction_Matlab - Calls the function that has been set with 4441 TSSetFunctionMatlab(). 4442 4443 Collective on TS 4444 4445 Input Parameters: 4446 + snes - the TS context 4447 - u - input vector 4448 4449 Output Parameter: 4450 . y - function vector, as set by TSSetFunction() 4451 4452 Notes: 4453 TSComputeFunction() is typically used within nonlinear solvers 4454 implementations, so most users would not generally call this routine 4455 themselves. 4456 4457 Level: developer 4458 4459 .keywords: TS, nonlinear, compute, function 4460 4461 .seealso: TSSetFunction(), TSGetFunction() 4462 */ 4463 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx) 4464 { 4465 PetscErrorCode ierr; 4466 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4467 int nlhs = 1,nrhs = 7; 4468 mxArray *plhs[1],*prhs[7]; 4469 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 4470 4471 PetscFunctionBegin; 4472 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 4473 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4474 PetscValidHeaderSpecific(udot,VEC_CLASSID,4); 4475 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 4476 PetscCheckSameComm(snes,1,u,3); 4477 PetscCheckSameComm(snes,1,y,5); 4478 4479 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4480 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4481 ierr = PetscMemcpy(&lxdot,&udot,sizeof(udot));CHKERRQ(ierr); 4482 ierr = PetscMemcpy(&ly,&y,sizeof(u));CHKERRQ(ierr); 4483 4484 prhs[0] = mxCreateDoubleScalar((double)ls); 4485 prhs[1] = mxCreateDoubleScalar(time); 4486 prhs[2] = mxCreateDoubleScalar((double)lx); 4487 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4488 prhs[4] = mxCreateDoubleScalar((double)ly); 4489 prhs[5] = mxCreateString(sctx->funcname); 4490 prhs[6] = sctx->ctx; 4491 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 4492 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4493 mxDestroyArray(prhs[0]); 4494 mxDestroyArray(prhs[1]); 4495 mxDestroyArray(prhs[2]); 4496 mxDestroyArray(prhs[3]); 4497 mxDestroyArray(prhs[4]); 4498 mxDestroyArray(prhs[5]); 4499 mxDestroyArray(plhs[0]); 4500 PetscFunctionReturn(0); 4501 } 4502 4503 4504 #undef __FUNCT__ 4505 #define __FUNCT__ "TSSetFunctionMatlab" 4506 /* 4507 TSSetFunctionMatlab - Sets the function evaluation routine and function 4508 vector for use by the TS routines in solving ODEs 4509 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4510 4511 Logically Collective on TS 4512 4513 Input Parameters: 4514 + ts - the TS context 4515 - func - function evaluation routine 4516 4517 Calling sequence of func: 4518 $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx); 4519 4520 Level: beginner 4521 4522 .keywords: TS, nonlinear, set, function 4523 4524 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4525 */ 4526 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 4527 { 4528 PetscErrorCode ierr; 4529 TSMatlabContext *sctx; 4530 4531 PetscFunctionBegin; 4532 /* currently sctx is memory bleed */ 4533 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4534 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4535 /* 4536 This should work, but it doesn't 4537 sctx->ctx = ctx; 4538 mexMakeArrayPersistent(sctx->ctx); 4539 */ 4540 sctx->ctx = mxDuplicateArray(ctx); 4541 4542 ierr = TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4543 PetscFunctionReturn(0); 4544 } 4545 4546 #undef __FUNCT__ 4547 #define __FUNCT__ "TSComputeJacobian_Matlab" 4548 /* 4549 TSComputeJacobian_Matlab - Calls the function that has been set with 4550 TSSetJacobianMatlab(). 4551 4552 Collective on TS 4553 4554 Input Parameters: 4555 + ts - the TS context 4556 . u - input vector 4557 . A, B - the matrices 4558 - ctx - user context 4559 4560 Output Parameter: 4561 . flag - structure of the matrix 4562 4563 Level: developer 4564 4565 .keywords: TS, nonlinear, compute, function 4566 4567 .seealso: TSSetFunction(), TSGetFunction() 4568 @*/ 4569 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4570 { 4571 PetscErrorCode ierr; 4572 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4573 int nlhs = 2,nrhs = 9; 4574 mxArray *plhs[2],*prhs[9]; 4575 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 4576 4577 PetscFunctionBegin; 4578 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4579 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4580 4581 /* call Matlab function in ctx with arguments u and y */ 4582 4583 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4584 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4585 ierr = PetscMemcpy(&lxdot,&udot,sizeof(u));CHKERRQ(ierr); 4586 ierr = PetscMemcpy(&lA,A,sizeof(u));CHKERRQ(ierr); 4587 ierr = PetscMemcpy(&lB,B,sizeof(u));CHKERRQ(ierr); 4588 4589 prhs[0] = mxCreateDoubleScalar((double)ls); 4590 prhs[1] = mxCreateDoubleScalar((double)time); 4591 prhs[2] = mxCreateDoubleScalar((double)lx); 4592 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4593 prhs[4] = mxCreateDoubleScalar((double)shift); 4594 prhs[5] = mxCreateDoubleScalar((double)lA); 4595 prhs[6] = mxCreateDoubleScalar((double)lB); 4596 prhs[7] = mxCreateString(sctx->funcname); 4597 prhs[8] = sctx->ctx; 4598 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 4599 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4600 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4601 mxDestroyArray(prhs[0]); 4602 mxDestroyArray(prhs[1]); 4603 mxDestroyArray(prhs[2]); 4604 mxDestroyArray(prhs[3]); 4605 mxDestroyArray(prhs[4]); 4606 mxDestroyArray(prhs[5]); 4607 mxDestroyArray(prhs[6]); 4608 mxDestroyArray(prhs[7]); 4609 mxDestroyArray(plhs[0]); 4610 mxDestroyArray(plhs[1]); 4611 PetscFunctionReturn(0); 4612 } 4613 4614 4615 #undef __FUNCT__ 4616 #define __FUNCT__ "TSSetJacobianMatlab" 4617 /* 4618 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4619 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 4620 4621 Logically Collective on TS 4622 4623 Input Parameters: 4624 + ts - the TS context 4625 . A,B - Jacobian matrices 4626 . func - function evaluation routine 4627 - ctx - user context 4628 4629 Calling sequence of func: 4630 $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx); 4631 4632 4633 Level: developer 4634 4635 .keywords: TS, nonlinear, set, function 4636 4637 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4638 */ 4639 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 4640 { 4641 PetscErrorCode ierr; 4642 TSMatlabContext *sctx; 4643 4644 PetscFunctionBegin; 4645 /* currently sctx is memory bleed */ 4646 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4647 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4648 /* 4649 This should work, but it doesn't 4650 sctx->ctx = ctx; 4651 mexMakeArrayPersistent(sctx->ctx); 4652 */ 4653 sctx->ctx = mxDuplicateArray(ctx); 4654 4655 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4656 PetscFunctionReturn(0); 4657 } 4658 4659 #undef __FUNCT__ 4660 #define __FUNCT__ "TSMonitor_Matlab" 4661 /* 4662 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 4663 4664 Collective on TS 4665 4666 .seealso: TSSetFunction(), TSGetFunction() 4667 @*/ 4668 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx) 4669 { 4670 PetscErrorCode ierr; 4671 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4672 int nlhs = 1,nrhs = 6; 4673 mxArray *plhs[1],*prhs[6]; 4674 long long int lx = 0,ls = 0; 4675 4676 PetscFunctionBegin; 4677 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4678 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 4679 4680 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4681 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4682 4683 prhs[0] = mxCreateDoubleScalar((double)ls); 4684 prhs[1] = mxCreateDoubleScalar((double)it); 4685 prhs[2] = mxCreateDoubleScalar((double)time); 4686 prhs[3] = mxCreateDoubleScalar((double)lx); 4687 prhs[4] = mxCreateString(sctx->funcname); 4688 prhs[5] = sctx->ctx; 4689 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 4690 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4691 mxDestroyArray(prhs[0]); 4692 mxDestroyArray(prhs[1]); 4693 mxDestroyArray(prhs[2]); 4694 mxDestroyArray(prhs[3]); 4695 mxDestroyArray(prhs[4]); 4696 mxDestroyArray(plhs[0]); 4697 PetscFunctionReturn(0); 4698 } 4699 4700 4701 #undef __FUNCT__ 4702 #define __FUNCT__ "TSMonitorSetMatlab" 4703 /* 4704 TSMonitorSetMatlab - Sets the monitor function from Matlab 4705 4706 Level: developer 4707 4708 .keywords: TS, nonlinear, set, function 4709 4710 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4711 */ 4712 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 4713 { 4714 PetscErrorCode ierr; 4715 TSMatlabContext *sctx; 4716 4717 PetscFunctionBegin; 4718 /* currently sctx is memory bleed */ 4719 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4720 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4721 /* 4722 This should work, but it doesn't 4723 sctx->ctx = ctx; 4724 mexMakeArrayPersistent(sctx->ctx); 4725 */ 4726 sctx->ctx = mxDuplicateArray(ctx); 4727 4728 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 4729 PetscFunctionReturn(0); 4730 } 4731 #endif 4732 4733 4734 4735 #undef __FUNCT__ 4736 #define __FUNCT__ "TSMonitorLGSolution" 4737 /*@C 4738 TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector 4739 in a time based line graph 4740 4741 Collective on TS 4742 4743 Input Parameters: 4744 + ts - the TS context 4745 . step - current time-step 4746 . ptime - current time 4747 - lg - a line graph object 4748 4749 Level: intermediate 4750 4751 Notes: each process in a parallel run displays its component solutions in a separate window 4752 4753 .keywords: TS, vector, monitor, view 4754 4755 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4756 @*/ 4757 PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4758 { 4759 PetscErrorCode ierr; 4760 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4761 const PetscScalar *yy; 4762 PetscInt dim; 4763 4764 PetscFunctionBegin; 4765 if (!step) { 4766 PetscDrawAxis axis; 4767 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4768 ierr = PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");CHKERRQ(ierr); 4769 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4770 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4771 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4772 } 4773 ierr = VecGetArrayRead(u,&yy);CHKERRQ(ierr); 4774 #if defined(PETSC_USE_COMPLEX) 4775 { 4776 PetscReal *yreal; 4777 PetscInt i,n; 4778 ierr = VecGetLocalSize(u,&n);CHKERRQ(ierr); 4779 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 4780 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4781 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4782 ierr = PetscFree(yreal);CHKERRQ(ierr); 4783 } 4784 #else 4785 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4786 #endif 4787 ierr = VecRestoreArrayRead(u,&yy);CHKERRQ(ierr); 4788 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 4789 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4790 } 4791 PetscFunctionReturn(0); 4792 } 4793 4794 #undef __FUNCT__ 4795 #define __FUNCT__ "TSMonitorLGError" 4796 /*@C 4797 TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector 4798 in a time based line graph 4799 4800 Collective on TS 4801 4802 Input Parameters: 4803 + ts - the TS context 4804 . step - current time-step 4805 . ptime - current time 4806 - lg - a line graph object 4807 4808 Level: intermediate 4809 4810 Notes: 4811 Only for sequential solves. 4812 4813 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 4814 4815 Options Database Keys: 4816 . -ts_monitor_lg_error - create a graphical monitor of error history 4817 4818 .keywords: TS, vector, monitor, view 4819 4820 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 4821 @*/ 4822 PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4823 { 4824 PetscErrorCode ierr; 4825 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4826 const PetscScalar *yy; 4827 Vec y; 4828 PetscInt dim; 4829 4830 PetscFunctionBegin; 4831 if (!step) { 4832 PetscDrawAxis axis; 4833 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4834 ierr = PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");CHKERRQ(ierr); 4835 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4836 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4837 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4838 } 4839 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 4840 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 4841 ierr = VecAXPY(y,-1.0,u);CHKERRQ(ierr); 4842 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 4843 #if defined(PETSC_USE_COMPLEX) 4844 { 4845 PetscReal *yreal; 4846 PetscInt i,n; 4847 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 4848 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 4849 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4850 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4851 ierr = PetscFree(yreal);CHKERRQ(ierr); 4852 } 4853 #else 4854 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4855 #endif 4856 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 4857 ierr = VecDestroy(&y);CHKERRQ(ierr); 4858 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 4859 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4860 } 4861 PetscFunctionReturn(0); 4862 } 4863 4864 #undef __FUNCT__ 4865 #define __FUNCT__ "TSMonitorLGSNESIterations" 4866 PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4867 { 4868 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4869 PetscReal x = ptime,y; 4870 PetscErrorCode ierr; 4871 PetscInt its; 4872 4873 PetscFunctionBegin; 4874 if (!n) { 4875 PetscDrawAxis axis; 4876 4877 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4878 ierr = PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");CHKERRQ(ierr); 4879 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4880 4881 ctx->snes_its = 0; 4882 } 4883 ierr = TSGetSNESIterations(ts,&its);CHKERRQ(ierr); 4884 y = its - ctx->snes_its; 4885 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4886 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 4887 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4888 } 4889 ctx->snes_its = its; 4890 PetscFunctionReturn(0); 4891 } 4892 4893 #undef __FUNCT__ 4894 #define __FUNCT__ "TSMonitorLGKSPIterations" 4895 PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4896 { 4897 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4898 PetscReal x = ptime,y; 4899 PetscErrorCode ierr; 4900 PetscInt its; 4901 4902 PetscFunctionBegin; 4903 if (!n) { 4904 PetscDrawAxis axis; 4905 4906 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4907 ierr = PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");CHKERRQ(ierr); 4908 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4909 4910 ctx->ksp_its = 0; 4911 } 4912 ierr = TSGetKSPIterations(ts,&its);CHKERRQ(ierr); 4913 y = its - ctx->ksp_its; 4914 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4915 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 4916 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4917 } 4918 ctx->ksp_its = its; 4919 PetscFunctionReturn(0); 4920 } 4921 4922 #undef __FUNCT__ 4923 #define __FUNCT__ "TSComputeLinearStability" 4924 /*@ 4925 TSComputeLinearStability - computes the linear stability function at a point 4926 4927 Collective on TS and Vec 4928 4929 Input Parameters: 4930 + ts - the TS context 4931 - xr,xi - real and imaginary part of input arguments 4932 4933 Output Parameters: 4934 . yr,yi - real and imaginary part of function value 4935 4936 Level: developer 4937 4938 .keywords: TS, compute 4939 4940 .seealso: TSSetRHSFunction(), TSComputeIFunction() 4941 @*/ 4942 PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 4943 { 4944 PetscErrorCode ierr; 4945 4946 PetscFunctionBegin; 4947 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4948 if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method"); 4949 ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr); 4950 PetscFunctionReturn(0); 4951 } 4952