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