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