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