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