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