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_save_trajectories - checkpoint the solution at each time-step 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 <xleft,yleft,xright,yright> - Monitor solution graphically with phase diagram, requires problem with exactly 2 degrees of freedom 55 . -ts_monitor_draw_error - Monitor error graphically, requires use to have provided TSSetSolutionFunction() 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_trajectory","Save the solution at each timestep","TSSetSaveTrajectory",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 ierr = TSTrajectorySetType(ts->trajectory,TSTRAJECTORYBASIC);CHKERRQ(ierr); 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__ "TSAdjointGetCostGradients" 1696 /*@ 1697 TSAdjointGetCostGradients - 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 + lambda - vectors containing the gradients of the cost functions with respect to the ODE/DAE solution variables 1706 - mu - vectors containing the gradients of the cost functions with respect to the problem parameters 1707 1708 Level: intermediate 1709 1710 .seealso: TSGetTimeStep() 1711 1712 .keywords: TS, timestep, get, sensitivity 1713 @*/ 1714 PetscErrorCode TSAdjointGetCostGradients(TS ts,PetscInt *numcost,Vec **lambda,Vec **mu) 1715 { 1716 PetscFunctionBegin; 1717 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1718 if (numcost) *numcost = ts->numcost; 1719 if (lambda) *lambda = ts->vecs_sensi; 1720 if (mu) *mu = 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 TSAdjointSetCostGradients() first"); 1915 1916 ierr = VecDuplicateVecs(ts->vecs_sensi[0],ts->numcost,&ts->vecs_drdy);CHKERRQ(ierr); 1917 if (ts->vecs_sensip){ 1918 ierr = VecDuplicateVecs(ts->vecs_sensip[0],ts->numcost,&ts->vecs_drdp);CHKERRQ(ierr); 1919 } 1920 1921 if (ts->ops->adjointsetup) { 1922 ierr = (*ts->ops->adjointsetup)(ts);CHKERRQ(ierr); 1923 } 1924 ts->adjointsetupcalled = PETSC_TRUE; 1925 PetscFunctionReturn(0); 1926 } 1927 1928 #undef __FUNCT__ 1929 #define __FUNCT__ "TSReset" 1930 /*@ 1931 TSReset - Resets a TS context and removes any allocated Vecs and Mats. 1932 1933 Collective on TS 1934 1935 Input Parameter: 1936 . ts - the TS context obtained from TSCreate() 1937 1938 Level: beginner 1939 1940 .keywords: TS, timestep, reset 1941 1942 .seealso: TSCreate(), TSSetup(), TSDestroy() 1943 @*/ 1944 PetscErrorCode TSReset(TS ts) 1945 { 1946 PetscErrorCode ierr; 1947 1948 PetscFunctionBegin; 1949 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1950 1951 if (ts->ops->reset) { 1952 ierr = (*ts->ops->reset)(ts);CHKERRQ(ierr); 1953 } 1954 if (ts->snes) {ierr = SNESReset(ts->snes);CHKERRQ(ierr);} 1955 1956 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 1957 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 1958 ierr = VecDestroy(&ts->Frhs);CHKERRQ(ierr); 1959 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1960 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 1961 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 1962 ierr = VecDestroyVecs(ts->nwork,&ts->work);CHKERRQ(ierr); 1963 ierr = VecDestroyVecs(ts->numcost,&ts->vecs_drdy);CHKERRQ(ierr); 1964 if (ts->vecs_drdp){ 1965 ierr = VecDestroyVecs(ts->numcost,&ts->vecs_drdp);CHKERRQ(ierr); 1966 } 1967 ts->vecs_sensi = NULL; 1968 ts->vecs_sensip = NULL; 1969 ierr = MatDestroy(&ts->Jacp);CHKERRQ(ierr); 1970 ierr = VecDestroy(&ts->vec_costintegral);CHKERRQ(ierr); 1971 ierr = VecDestroy(&ts->vec_costintegrand);CHKERRQ(ierr); 1972 ts->setupcalled = PETSC_FALSE; 1973 PetscFunctionReturn(0); 1974 } 1975 1976 #undef __FUNCT__ 1977 #define __FUNCT__ "TSDestroy" 1978 /*@ 1979 TSDestroy - Destroys the timestepper context that was created 1980 with TSCreate(). 1981 1982 Collective on TS 1983 1984 Input Parameter: 1985 . ts - the TS context obtained from TSCreate() 1986 1987 Level: beginner 1988 1989 .keywords: TS, timestepper, destroy 1990 1991 .seealso: TSCreate(), TSSetUp(), TSSolve() 1992 @*/ 1993 PetscErrorCode TSDestroy(TS *ts) 1994 { 1995 PetscErrorCode ierr; 1996 1997 PetscFunctionBegin; 1998 if (!*ts) PetscFunctionReturn(0); 1999 PetscValidHeaderSpecific((*ts),TS_CLASSID,1); 2000 if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; PetscFunctionReturn(0);} 2001 2002 ierr = TSReset((*ts));CHKERRQ(ierr); 2003 2004 /* if memory was published with SAWs then destroy it */ 2005 ierr = PetscObjectSAWsViewOff((PetscObject)*ts);CHKERRQ(ierr); 2006 if ((*ts)->ops->destroy) {ierr = (*(*ts)->ops->destroy)((*ts));CHKERRQ(ierr);} 2007 2008 ierr = TSTrajectoryDestroy(&(*ts)->trajectory);CHKERRQ(ierr); 2009 2010 ierr = TSAdaptDestroy(&(*ts)->adapt);CHKERRQ(ierr); 2011 if ((*ts)->event) { 2012 ierr = TSEventMonitorDestroy(&(*ts)->event);CHKERRQ(ierr); 2013 } 2014 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 2015 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 2016 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 2017 2018 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 2019 PetscFunctionReturn(0); 2020 } 2021 2022 #undef __FUNCT__ 2023 #define __FUNCT__ "TSGetSNES" 2024 /*@ 2025 TSGetSNES - Returns the SNES (nonlinear solver) associated with 2026 a TS (timestepper) context. Valid only for nonlinear problems. 2027 2028 Not Collective, but SNES is parallel if TS is parallel 2029 2030 Input Parameter: 2031 . ts - the TS context obtained from TSCreate() 2032 2033 Output Parameter: 2034 . snes - the nonlinear solver context 2035 2036 Notes: 2037 The user can then directly manipulate the SNES context to set various 2038 options, etc. Likewise, the user can then extract and manipulate the 2039 KSP, KSP, and PC contexts as well. 2040 2041 TSGetSNES() does not work for integrators that do not use SNES; in 2042 this case TSGetSNES() returns NULL in snes. 2043 2044 Level: beginner 2045 2046 .keywords: timestep, get, SNES 2047 @*/ 2048 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 2049 { 2050 PetscErrorCode ierr; 2051 2052 PetscFunctionBegin; 2053 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2054 PetscValidPointer(snes,2); 2055 if (!ts->snes) { 2056 ierr = SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);CHKERRQ(ierr); 2057 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 2058 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->snes);CHKERRQ(ierr); 2059 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 2060 if (ts->dm) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 2061 if (ts->problem_type == TS_LINEAR) { 2062 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 2063 } 2064 } 2065 *snes = ts->snes; 2066 PetscFunctionReturn(0); 2067 } 2068 2069 #undef __FUNCT__ 2070 #define __FUNCT__ "TSSetSNES" 2071 /*@ 2072 TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context 2073 2074 Collective 2075 2076 Input Parameter: 2077 + ts - the TS context obtained from TSCreate() 2078 - snes - the nonlinear solver context 2079 2080 Notes: 2081 Most users should have the TS created by calling TSGetSNES() 2082 2083 Level: developer 2084 2085 .keywords: timestep, set, SNES 2086 @*/ 2087 PetscErrorCode TSSetSNES(TS ts,SNES snes) 2088 { 2089 PetscErrorCode ierr; 2090 PetscErrorCode (*func)(SNES,Vec,Mat,Mat,void*); 2091 2092 PetscFunctionBegin; 2093 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2094 PetscValidHeaderSpecific(snes,SNES_CLASSID,2); 2095 ierr = PetscObjectReference((PetscObject)snes);CHKERRQ(ierr); 2096 ierr = SNESDestroy(&ts->snes);CHKERRQ(ierr); 2097 2098 ts->snes = snes; 2099 2100 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 2101 ierr = SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);CHKERRQ(ierr); 2102 if (func == SNESTSFormJacobian) { 2103 ierr = SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);CHKERRQ(ierr); 2104 } 2105 PetscFunctionReturn(0); 2106 } 2107 2108 #undef __FUNCT__ 2109 #define __FUNCT__ "TSGetKSP" 2110 /*@ 2111 TSGetKSP - Returns the KSP (linear solver) associated with 2112 a TS (timestepper) context. 2113 2114 Not Collective, but KSP is parallel if TS is parallel 2115 2116 Input Parameter: 2117 . ts - the TS context obtained from TSCreate() 2118 2119 Output Parameter: 2120 . ksp - the nonlinear solver context 2121 2122 Notes: 2123 The user can then directly manipulate the KSP context to set various 2124 options, etc. Likewise, the user can then extract and manipulate the 2125 KSP and PC contexts as well. 2126 2127 TSGetKSP() does not work for integrators that do not use KSP; 2128 in this case TSGetKSP() returns NULL in ksp. 2129 2130 Level: beginner 2131 2132 .keywords: timestep, get, KSP 2133 @*/ 2134 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 2135 { 2136 PetscErrorCode ierr; 2137 SNES snes; 2138 2139 PetscFunctionBegin; 2140 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2141 PetscValidPointer(ksp,2); 2142 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 2143 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 2144 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2145 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 2146 PetscFunctionReturn(0); 2147 } 2148 2149 /* ----------- Routines to set solver parameters ---------- */ 2150 2151 #undef __FUNCT__ 2152 #define __FUNCT__ "TSGetDuration" 2153 /*@ 2154 TSGetDuration - Gets the maximum number of timesteps to use and 2155 maximum time for iteration. 2156 2157 Not Collective 2158 2159 Input Parameters: 2160 + ts - the TS context obtained from TSCreate() 2161 . maxsteps - maximum number of iterations to use, or NULL 2162 - maxtime - final time to iterate to, or NULL 2163 2164 Level: intermediate 2165 2166 .keywords: TS, timestep, get, maximum, iterations, time 2167 @*/ 2168 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 2169 { 2170 PetscFunctionBegin; 2171 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2172 if (maxsteps) { 2173 PetscValidIntPointer(maxsteps,2); 2174 *maxsteps = ts->max_steps; 2175 } 2176 if (maxtime) { 2177 PetscValidScalarPointer(maxtime,3); 2178 *maxtime = ts->max_time; 2179 } 2180 PetscFunctionReturn(0); 2181 } 2182 2183 #undef __FUNCT__ 2184 #define __FUNCT__ "TSSetDuration" 2185 /*@ 2186 TSSetDuration - Sets the maximum number of timesteps to use and 2187 maximum time for iteration. 2188 2189 Logically Collective on TS 2190 2191 Input Parameters: 2192 + ts - the TS context obtained from TSCreate() 2193 . maxsteps - maximum number of iterations to use 2194 - maxtime - final time to iterate to 2195 2196 Options Database Keys: 2197 . -ts_max_steps <maxsteps> - Sets maxsteps 2198 . -ts_final_time <maxtime> - Sets maxtime 2199 2200 Notes: 2201 The default maximum number of iterations is 5000. Default time is 5.0 2202 2203 Level: intermediate 2204 2205 .keywords: TS, timestep, set, maximum, iterations 2206 2207 .seealso: TSSetExactFinalTime() 2208 @*/ 2209 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 2210 { 2211 PetscFunctionBegin; 2212 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2213 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 2214 PetscValidLogicalCollectiveReal(ts,maxtime,2); 2215 if (maxsteps >= 0) ts->max_steps = maxsteps; 2216 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 2217 PetscFunctionReturn(0); 2218 } 2219 2220 #undef __FUNCT__ 2221 #define __FUNCT__ "TSSetSolution" 2222 /*@ 2223 TSSetSolution - Sets the initial solution vector 2224 for use by the TS routines. 2225 2226 Logically Collective on TS and Vec 2227 2228 Input Parameters: 2229 + ts - the TS context obtained from TSCreate() 2230 - u - the solution vector 2231 2232 Level: beginner 2233 2234 .keywords: TS, timestep, set, solution, initial conditions 2235 @*/ 2236 PetscErrorCode TSSetSolution(TS ts,Vec u) 2237 { 2238 PetscErrorCode ierr; 2239 DM dm; 2240 2241 PetscFunctionBegin; 2242 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2243 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2244 ierr = PetscObjectReference((PetscObject)u);CHKERRQ(ierr); 2245 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 2246 2247 ts->vec_sol = u; 2248 2249 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2250 ierr = DMShellSetGlobalVector(dm,u);CHKERRQ(ierr); 2251 PetscFunctionReturn(0); 2252 } 2253 2254 #undef __FUNCT__ 2255 #define __FUNCT__ "TSAdjointSetSteps" 2256 /*@ 2257 TSAdjointSetSteps - Sets the number of steps the adjoint solver should take backward in time 2258 2259 Logically Collective on TS 2260 2261 Input Parameters: 2262 + ts - the TS context obtained from TSCreate() 2263 . steps - number of steps to use 2264 2265 Level: intermediate 2266 2267 Notes: Normally one does not call this and TSAdjointSolve() integrates back to the original timestep. One can call this 2268 so as to integrate back to less than the original timestep 2269 2270 .keywords: TS, timestep, set, maximum, iterations 2271 2272 .seealso: TSSetExactFinalTime() 2273 @*/ 2274 PetscErrorCode TSAdjointSetSteps(TS ts,PetscInt steps) 2275 { 2276 PetscFunctionBegin; 2277 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2278 PetscValidLogicalCollectiveInt(ts,steps,2); 2279 if (steps < 0) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back a negative number of steps"); 2280 if (steps > ts->total_steps) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"Cannot step back more than the total number of forward steps"); 2281 ts->adjoint_max_steps = steps; 2282 PetscFunctionReturn(0); 2283 } 2284 2285 #undef __FUNCT__ 2286 #define __FUNCT__ "TSAdjointSetCostGradients" 2287 /*@ 2288 TSAdjointSetCostGradients - Sets the initial value of the gradients of the cost function w.r.t. initial conditions and w.r.t. the problem parameters 2289 for use by the TSAdjoint routines. 2290 2291 Logically Collective on TS and Vec 2292 2293 Input Parameters: 2294 + ts - the TS context obtained from TSCreate() 2295 . lambda - gradients with respect to the initial condition variables, the dimension and parallel layout of these vectors is the same as the ODE solution vector 2296 - mu - gradients with respect to the parameters, the number of entries in these vectors is the same as the number of parameters 2297 2298 Level: beginner 2299 2300 Notes: the entries in these vectors must be correctly initialized with the values lamda_i = df/dy|finaltime mu_i = df/dp|finaltime 2301 2302 .keywords: TS, timestep, set, sensitivity, initial conditions 2303 @*/ 2304 PetscErrorCode TSAdjointSetCostGradients(TS ts,PetscInt numcost,Vec *lambda,Vec *mu) 2305 { 2306 PetscFunctionBegin; 2307 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2308 PetscValidPointer(lambda,2); 2309 ts->vecs_sensi = lambda; 2310 ts->vecs_sensip = mu; 2311 ts->numcost = numcost; 2312 PetscFunctionReturn(0); 2313 } 2314 2315 #undef __FUNCT__ 2316 #define __FUNCT__ "TSAdjointSetRHSJacobian" 2317 /*@C 2318 TSAdjointSetRHSJacobian - Sets the function that computes the Jacobian of G w.r.t. the parameters p where y_t = G(y,p,t), as well as the location to store the matrix. 2319 2320 Logically Collective on TS 2321 2322 Input Parameters: 2323 + ts - The TS context obtained from TSCreate() 2324 - func - The function 2325 2326 Calling sequence of func: 2327 $ func (TS ts,PetscReal t,Vec y,Mat A,void *ctx); 2328 + t - current timestep 2329 . y - input vector (current ODE solution) 2330 . A - output matrix 2331 - ctx - [optional] user-defined function context 2332 2333 Level: intermediate 2334 2335 Notes: Amat has the same number of rows and the same row parallel layout as u, Amat has the same number of columns and parallel layout as p 2336 2337 .keywords: TS, sensitivity 2338 .seealso: 2339 @*/ 2340 PetscErrorCode TSAdjointSetRHSJacobian(TS ts,Mat Amat,PetscErrorCode (*func)(TS,PetscReal,Vec,Mat,void*),void *ctx) 2341 { 2342 PetscErrorCode ierr; 2343 2344 PetscFunctionBegin; 2345 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2346 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 2347 2348 ts->rhsjacobianp = func; 2349 ts->rhsjacobianpctx = ctx; 2350 if(Amat) { 2351 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 2352 ierr = MatDestroy(&ts->Jacp);CHKERRQ(ierr); 2353 ts->Jacp = Amat; 2354 } 2355 PetscFunctionReturn(0); 2356 } 2357 2358 #undef __FUNCT__ 2359 #define __FUNCT__ "TSAdjointComputeRHSJacobian" 2360 /*@C 2361 TSAdjointComputeRHSJacobian - Runs the user-defined Jacobian function. 2362 2363 Collective on TS 2364 2365 Input Parameters: 2366 . ts - The TS context obtained from TSCreate() 2367 2368 Level: developer 2369 2370 .keywords: TS, sensitivity 2371 .seealso: TSAdjointSetRHSJacobian() 2372 @*/ 2373 PetscErrorCode TSAdjointComputeRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat) 2374 { 2375 PetscErrorCode ierr; 2376 2377 PetscFunctionBegin; 2378 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2379 PetscValidHeaderSpecific(X,VEC_CLASSID,3); 2380 PetscValidPointer(Amat,4); 2381 2382 PetscStackPush("TS user JacobianP function for sensitivity analysis"); 2383 ierr = (*ts->rhsjacobianp)(ts,t,X,Amat,ts->rhsjacobianpctx); CHKERRQ(ierr); 2384 PetscStackPop; 2385 PetscFunctionReturn(0); 2386 } 2387 2388 #undef __FUNCT__ 2389 #define __FUNCT__ "TSAdjointSetCostIntegrand" 2390 /*@C 2391 TSAdjointSetCostIntegrand - Sets the routine for evaluating the integral term in one or more cost functions 2392 2393 Logically Collective on TS 2394 2395 Input Parameters: 2396 + ts - the TS context obtained from TSCreate() 2397 . numcost - number of gradients to be computed, this is the number of cost functions 2398 . rf - routine for evaluating the integrand function 2399 . drdyf - function that computes the gradients of the r's with respect to y,NULL if not a function y 2400 . drdpf - function that computes the gradients of the r's with respect to p, NULL if not a function of p 2401 - ctx - [optional] user-defined context for private data for the function evaluation routine (may be NULL) 2402 2403 Calling sequence of rf: 2404 $ rf(TS ts,PetscReal t,Vec y,Vec f[],void *ctx); 2405 2406 + t - current timestep 2407 . y - input vector 2408 . f - function result; one vector entry for each cost function 2409 - ctx - [optional] user-defined function context 2410 2411 Calling sequence of drdyf: 2412 $ PetscErroCode drdyf(TS ts,PetscReal t,Vec y,Vec *drdy,void *ctx); 2413 2414 Calling sequence of drdpf: 2415 $ PetscErroCode drdpf(TS ts,PetscReal t,Vec y,Vec *drdp,void *ctx); 2416 2417 Level: intermediate 2418 2419 Notes: For optimization there is generally a single cost function, numcost = 1. For sensitivities there may be multiple cost functions 2420 2421 .keywords: TS, sensitivity analysis, timestep, set, quadrature, function 2422 2423 .seealso: TSAdjointSetRHSJacobian(),TSAdjointGetCostGradients(), TSAdjointSetCostGradients() 2424 @*/ 2425 PetscErrorCode TSAdjointSetCostIntegrand(TS ts,PetscInt numcost, PetscErrorCode (*rf)(TS,PetscReal,Vec,Vec,void*), 2426 PetscErrorCode (*drdyf)(TS,PetscReal,Vec,Vec*,void*), 2427 PetscErrorCode (*drdpf)(TS,PetscReal,Vec,Vec*,void*),void *ctx) 2428 { 2429 PetscErrorCode ierr; 2430 2431 PetscFunctionBegin; 2432 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2433 if (!ts->numcost) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Call TSAdjointSetCostGradients() first so that the number of cost functions can be determined."); 2434 if (ts->numcost && ts->numcost!=numcost) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"The number of cost functions (2rd parameter of TSAdjointSetCostIntegrand()) is inconsistent with the one set by TSAdjointSetCostGradients()"); 2435 2436 ierr = VecCreateSeq(PETSC_COMM_SELF,numcost,&ts->vec_costintegral);CHKERRQ(ierr); 2437 ierr = VecDuplicate(ts->vec_costintegral,&ts->vec_costintegrand);CHKERRQ(ierr); 2438 ts->costintegrand = rf; 2439 ts->costintegrandctx = ctx; 2440 ts->drdyfunction = drdyf; 2441 ts->drdpfunction = drdpf; 2442 PetscFunctionReturn(0); 2443 } 2444 2445 #undef __FUNCT__ 2446 #define __FUNCT__ "TSAdjointGetCostIntegral" 2447 /*@ 2448 TSAdjointGetCostIntegral - Returns the values of the integral term in the cost functions. 2449 It is valid to call the routine after a backward run. 2450 2451 Not Collective 2452 2453 Input Parameter: 2454 . ts - the TS context obtained from TSCreate() 2455 2456 Output Parameter: 2457 . v - the vector containing the integrals for each cost function 2458 2459 Level: intermediate 2460 2461 .seealso: TSAdjointSetCostIntegrand() 2462 2463 .keywords: TS, sensitivity analysis 2464 @*/ 2465 PetscErrorCode TSAdjointGetCostIntegral(TS ts,Vec *v) 2466 { 2467 PetscFunctionBegin; 2468 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2469 PetscValidPointer(v,2); 2470 *v = ts->vec_costintegral; 2471 PetscFunctionReturn(0); 2472 } 2473 2474 #undef __FUNCT__ 2475 #define __FUNCT__ "TSAdjointComputeCostIntegrand" 2476 /*@ 2477 TSAdjointComputeCostIntegrand - Evaluates the integral function in the cost functions. 2478 2479 Input Parameters: 2480 + ts - the TS context 2481 . t - current time 2482 - y - state vector, i.e. current solution 2483 2484 Output Parameter: 2485 . q - vector of size numcost to hold the outputs 2486 2487 Note: 2488 Most users should not need to explicitly call this routine, as it 2489 is used internally within the sensitivity analysis context. 2490 2491 Level: developer 2492 2493 .keywords: TS, compute 2494 2495 .seealso: TSAdjointSetCostIntegrand() 2496 @*/ 2497 PetscErrorCode TSAdjointComputeCostIntegrand(TS ts,PetscReal t,Vec y,Vec q) 2498 { 2499 PetscErrorCode ierr; 2500 2501 PetscFunctionBegin; 2502 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2503 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2504 PetscValidHeaderSpecific(q,VEC_CLASSID,4); 2505 2506 ierr = PetscLogEventBegin(TS_FunctionEval,ts,y,q,0);CHKERRQ(ierr); 2507 if (ts->costintegrand) { 2508 PetscStackPush("TS user integrand in the cost function"); 2509 ierr = (*ts->costintegrand)(ts,t,y,q,ts->costintegrandctx);CHKERRQ(ierr); 2510 PetscStackPop; 2511 } else { 2512 ierr = VecZeroEntries(q);CHKERRQ(ierr); 2513 } 2514 2515 ierr = PetscLogEventEnd(TS_FunctionEval,ts,y,q,0);CHKERRQ(ierr); 2516 PetscFunctionReturn(0); 2517 } 2518 2519 #undef __FUNCT__ 2520 #define __FUNCT__ "TSAdjointComputeDRDYFunction" 2521 /*@ 2522 TSAdjointComputeDRDYFunction - Runs the user-defined DRDY function. 2523 2524 Collective on TS 2525 2526 Input Parameters: 2527 . ts - The TS context obtained from TSCreate() 2528 2529 Notes: 2530 TSAdjointComputeDRDYFunction() is typically used for sensitivity implementation, 2531 so most users would not generally call this routine themselves. 2532 2533 Level: developer 2534 2535 .keywords: TS, sensitivity 2536 .seealso: TSAdjointComputeDRDYFunction() 2537 @*/ 2538 PetscErrorCode TSAdjointComputeDRDYFunction(TS ts,PetscReal t,Vec y,Vec *drdy) 2539 { 2540 PetscErrorCode ierr; 2541 2542 PetscFunctionBegin; 2543 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2544 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2545 2546 PetscStackPush("TS user DRDY function for sensitivity analysis"); 2547 ierr = (*ts->drdyfunction)(ts,t,y,drdy,ts->costintegrandctx); CHKERRQ(ierr); 2548 PetscStackPop; 2549 PetscFunctionReturn(0); 2550 } 2551 2552 #undef __FUNCT__ 2553 #define __FUNCT__ "TSAdjointComputeDRDPFunction" 2554 /*@ 2555 TSAdjointComputeDRDPFunction - Runs the user-defined DRDP function. 2556 2557 Collective on TS 2558 2559 Input Parameters: 2560 . ts - The TS context obtained from TSCreate() 2561 2562 Notes: 2563 TSDRDPFunction() is typically used for sensitivity implementation, 2564 so most users would not generally call this routine themselves. 2565 2566 Level: developer 2567 2568 .keywords: TS, sensitivity 2569 .seealso: TSAdjointSetDRDPFunction() 2570 @*/ 2571 PetscErrorCode TSAdjointComputeDRDPFunction(TS ts,PetscReal t,Vec y,Vec *drdp) 2572 { 2573 PetscErrorCode ierr; 2574 2575 PetscFunctionBegin; 2576 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2577 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2578 2579 PetscStackPush("TS user DRDP function for sensitivity analysis"); 2580 ierr = (*ts->drdpfunction)(ts,t,y,drdp,ts->costintegrandctx); CHKERRQ(ierr); 2581 PetscStackPop; 2582 PetscFunctionReturn(0); 2583 } 2584 2585 #undef __FUNCT__ 2586 #define __FUNCT__ "TSSetPreStep" 2587 /*@C 2588 TSSetPreStep - Sets the general-purpose function 2589 called once at the beginning of each time step. 2590 2591 Logically Collective on TS 2592 2593 Input Parameters: 2594 + ts - The TS context obtained from TSCreate() 2595 - func - The function 2596 2597 Calling sequence of func: 2598 . func (TS ts); 2599 2600 Level: intermediate 2601 2602 Note: 2603 If a step is rejected, TSStep() will call this routine again before each attempt. 2604 The last completed time step number can be queried using TSGetTimeStepNumber(), the 2605 size of the step being attempted can be obtained using TSGetTimeStep(). 2606 2607 .keywords: TS, timestep 2608 .seealso: TSSetPreStage(), TSSetPostStage(), TSSetPostStep(), TSStep() 2609 @*/ 2610 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 2611 { 2612 PetscFunctionBegin; 2613 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2614 ts->prestep = func; 2615 PetscFunctionReturn(0); 2616 } 2617 2618 #undef __FUNCT__ 2619 #define __FUNCT__ "TSPreStep" 2620 /*@ 2621 TSPreStep - Runs the user-defined pre-step function. 2622 2623 Collective on TS 2624 2625 Input Parameters: 2626 . ts - The TS context obtained from TSCreate() 2627 2628 Notes: 2629 TSPreStep() is typically used within time stepping implementations, 2630 so most users would not generally call this routine themselves. 2631 2632 Level: developer 2633 2634 .keywords: TS, timestep 2635 .seealso: TSSetPreStep(), TSPreStage(), TSPostStage(), TSPostStep() 2636 @*/ 2637 PetscErrorCode TSPreStep(TS ts) 2638 { 2639 PetscErrorCode ierr; 2640 2641 PetscFunctionBegin; 2642 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2643 if (ts->prestep) { 2644 PetscStackCallStandard((*ts->prestep),(ts)); 2645 } 2646 PetscFunctionReturn(0); 2647 } 2648 2649 #undef __FUNCT__ 2650 #define __FUNCT__ "TSSetPreStage" 2651 /*@C 2652 TSSetPreStage - Sets the general-purpose function 2653 called once at the beginning of each stage. 2654 2655 Logically Collective on TS 2656 2657 Input Parameters: 2658 + ts - The TS context obtained from TSCreate() 2659 - func - The function 2660 2661 Calling sequence of func: 2662 . PetscErrorCode func(TS ts, PetscReal stagetime); 2663 2664 Level: intermediate 2665 2666 Note: 2667 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2668 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2669 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2670 2671 .keywords: TS, timestep 2672 .seealso: TSSetPostStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2673 @*/ 2674 PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal)) 2675 { 2676 PetscFunctionBegin; 2677 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2678 ts->prestage = func; 2679 PetscFunctionReturn(0); 2680 } 2681 2682 #undef __FUNCT__ 2683 #define __FUNCT__ "TSSetPostStage" 2684 /*@C 2685 TSSetPostStage - Sets the general-purpose function 2686 called once at the end of each stage. 2687 2688 Logically Collective on TS 2689 2690 Input Parameters: 2691 + ts - The TS context obtained from TSCreate() 2692 - func - The function 2693 2694 Calling sequence of func: 2695 . PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y); 2696 2697 Level: intermediate 2698 2699 Note: 2700 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2701 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2702 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2703 2704 .keywords: TS, timestep 2705 .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2706 @*/ 2707 PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS,PetscReal,PetscInt,Vec*)) 2708 { 2709 PetscFunctionBegin; 2710 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2711 ts->poststage = func; 2712 PetscFunctionReturn(0); 2713 } 2714 2715 #undef __FUNCT__ 2716 #define __FUNCT__ "TSPreStage" 2717 /*@ 2718 TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage() 2719 2720 Collective on TS 2721 2722 Input Parameters: 2723 . ts - The TS context obtained from TSCreate() 2724 stagetime - The absolute time of the current stage 2725 2726 Notes: 2727 TSPreStage() is typically used within time stepping implementations, 2728 most users would not generally call this routine themselves. 2729 2730 Level: developer 2731 2732 .keywords: TS, timestep 2733 .seealso: TSPostStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 2734 @*/ 2735 PetscErrorCode TSPreStage(TS ts, PetscReal stagetime) 2736 { 2737 PetscErrorCode ierr; 2738 2739 PetscFunctionBegin; 2740 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2741 if (ts->prestage) { 2742 PetscStackCallStandard((*ts->prestage),(ts,stagetime)); 2743 } 2744 PetscFunctionReturn(0); 2745 } 2746 2747 #undef __FUNCT__ 2748 #define __FUNCT__ "TSPostStage" 2749 /*@ 2750 TSPostStage - Runs the user-defined post-stage function set using TSSetPostStage() 2751 2752 Collective on TS 2753 2754 Input Parameters: 2755 . ts - The TS context obtained from TSCreate() 2756 stagetime - The absolute time of the current stage 2757 stageindex - Stage number 2758 Y - Array of vectors (of size = total number 2759 of stages) with the stage solutions 2760 2761 Notes: 2762 TSPostStage() is typically used within time stepping implementations, 2763 most users would not generally call this routine themselves. 2764 2765 Level: developer 2766 2767 .keywords: TS, timestep 2768 .seealso: TSPreStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 2769 @*/ 2770 PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y) 2771 { 2772 PetscErrorCode ierr; 2773 2774 PetscFunctionBegin; 2775 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2776 if (ts->poststage) { 2777 PetscStackCallStandard((*ts->poststage),(ts,stagetime,stageindex,Y)); 2778 } 2779 PetscFunctionReturn(0); 2780 } 2781 2782 #undef __FUNCT__ 2783 #define __FUNCT__ "TSSetPostStep" 2784 /*@C 2785 TSSetPostStep - Sets the general-purpose function 2786 called once at the end of each time step. 2787 2788 Logically Collective on TS 2789 2790 Input Parameters: 2791 + ts - The TS context obtained from TSCreate() 2792 - func - The function 2793 2794 Calling sequence of func: 2795 $ func (TS ts); 2796 2797 Level: intermediate 2798 2799 .keywords: TS, timestep 2800 .seealso: TSSetPreStep(), TSSetPreStage(), TSGetTimeStep(), TSGetTimeStepNumber(), TSGetTime() 2801 @*/ 2802 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 2803 { 2804 PetscFunctionBegin; 2805 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2806 ts->poststep = func; 2807 PetscFunctionReturn(0); 2808 } 2809 2810 #undef __FUNCT__ 2811 #define __FUNCT__ "TSPostStep" 2812 /*@ 2813 TSPostStep - Runs the user-defined post-step function. 2814 2815 Collective on TS 2816 2817 Input Parameters: 2818 . ts - The TS context obtained from TSCreate() 2819 2820 Notes: 2821 TSPostStep() is typically used within time stepping implementations, 2822 so most users would not generally call this routine themselves. 2823 2824 Level: developer 2825 2826 .keywords: TS, timestep 2827 @*/ 2828 PetscErrorCode TSPostStep(TS ts) 2829 { 2830 PetscErrorCode ierr; 2831 2832 PetscFunctionBegin; 2833 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2834 if (ts->poststep) { 2835 PetscStackCallStandard((*ts->poststep),(ts)); 2836 } 2837 PetscFunctionReturn(0); 2838 } 2839 2840 /* ------------ Routines to set performance monitoring options ----------- */ 2841 2842 #undef __FUNCT__ 2843 #define __FUNCT__ "TSMonitorSet" 2844 /*@C 2845 TSMonitorSet - Sets an ADDITIONAL function that is to be used at every 2846 timestep to display the iteration's progress. 2847 2848 Logically Collective on TS 2849 2850 Input Parameters: 2851 + ts - the TS context obtained from TSCreate() 2852 . monitor - monitoring routine 2853 . mctx - [optional] user-defined context for private data for the 2854 monitor routine (use NULL if no context is desired) 2855 - monitordestroy - [optional] routine that frees monitor context 2856 (may be NULL) 2857 2858 Calling sequence of monitor: 2859 $ int monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx) 2860 2861 + ts - the TS context 2862 . 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 2863 been interpolated to) 2864 . time - current time 2865 . u - current iterate 2866 - mctx - [optional] monitoring context 2867 2868 Notes: 2869 This routine adds an additional monitor to the list of monitors that 2870 already has been loaded. 2871 2872 Fortran notes: Only a single monitor function can be set for each TS object 2873 2874 Level: intermediate 2875 2876 .keywords: TS, timestep, set, monitor 2877 2878 .seealso: TSMonitorDefault(), TSMonitorCancel() 2879 @*/ 2880 PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**)) 2881 { 2882 PetscFunctionBegin; 2883 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2884 if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 2885 ts->monitor[ts->numbermonitors] = monitor; 2886 ts->monitordestroy[ts->numbermonitors] = mdestroy; 2887 ts->monitorcontext[ts->numbermonitors++] = (void*)mctx; 2888 PetscFunctionReturn(0); 2889 } 2890 2891 #undef __FUNCT__ 2892 #define __FUNCT__ "TSMonitorCancel" 2893 /*@C 2894 TSMonitorCancel - Clears all the monitors that have been set on a time-step object. 2895 2896 Logically Collective on TS 2897 2898 Input Parameters: 2899 . ts - the TS context obtained from TSCreate() 2900 2901 Notes: 2902 There is no way to remove a single, specific monitor. 2903 2904 Level: intermediate 2905 2906 .keywords: TS, timestep, set, monitor 2907 2908 .seealso: TSMonitorDefault(), TSMonitorSet() 2909 @*/ 2910 PetscErrorCode TSMonitorCancel(TS ts) 2911 { 2912 PetscErrorCode ierr; 2913 PetscInt i; 2914 2915 PetscFunctionBegin; 2916 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2917 for (i=0; i<ts->numbermonitors; i++) { 2918 if (ts->monitordestroy[i]) { 2919 ierr = (*ts->monitordestroy[i])(&ts->monitorcontext[i]);CHKERRQ(ierr); 2920 } 2921 } 2922 ts->numbermonitors = 0; 2923 PetscFunctionReturn(0); 2924 } 2925 2926 #undef __FUNCT__ 2927 #define __FUNCT__ "TSMonitorDefault" 2928 /*@ 2929 TSMonitorDefault - Sets the Default monitor 2930 2931 Level: intermediate 2932 2933 .keywords: TS, set, monitor 2934 2935 .seealso: TSMonitorDefault(), TSMonitorSet() 2936 @*/ 2937 PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy) 2938 { 2939 PetscErrorCode ierr; 2940 PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ts)); 2941 2942 PetscFunctionBegin; 2943 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2944 ierr = PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g\n",step,(double)ts->time_step,(double)ptime);CHKERRQ(ierr); 2945 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2946 PetscFunctionReturn(0); 2947 } 2948 2949 #undef __FUNCT__ 2950 #define __FUNCT__ "TSSetRetainStages" 2951 /*@ 2952 TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available. 2953 2954 Logically Collective on TS 2955 2956 Input Argument: 2957 . ts - time stepping context 2958 2959 Output Argument: 2960 . flg - PETSC_TRUE or PETSC_FALSE 2961 2962 Level: intermediate 2963 2964 .keywords: TS, set 2965 2966 .seealso: TSInterpolate(), TSSetPostStep() 2967 @*/ 2968 PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg) 2969 { 2970 PetscFunctionBegin; 2971 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2972 ts->retain_stages = flg; 2973 PetscFunctionReturn(0); 2974 } 2975 2976 #undef __FUNCT__ 2977 #define __FUNCT__ "TSInterpolate" 2978 /*@ 2979 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 2980 2981 Collective on TS 2982 2983 Input Argument: 2984 + ts - time stepping context 2985 - t - time to interpolate to 2986 2987 Output Argument: 2988 . U - state at given time 2989 2990 Notes: 2991 The user should call TSSetRetainStages() before taking a step in which interpolation will be requested. 2992 2993 Level: intermediate 2994 2995 Developer Notes: 2996 TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 2997 2998 .keywords: TS, set 2999 3000 .seealso: TSSetRetainStages(), TSSetPostStep() 3001 @*/ 3002 PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U) 3003 { 3004 PetscErrorCode ierr; 3005 3006 PetscFunctionBegin; 3007 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3008 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 3009 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); 3010 if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name); 3011 ierr = (*ts->ops->interpolate)(ts,t,U);CHKERRQ(ierr); 3012 PetscFunctionReturn(0); 3013 } 3014 3015 #undef __FUNCT__ 3016 #define __FUNCT__ "TSStep" 3017 /*@ 3018 TSStep - Steps one time step 3019 3020 Collective on TS 3021 3022 Input Parameter: 3023 . ts - the TS context obtained from TSCreate() 3024 3025 Level: developer 3026 3027 Notes: 3028 The public interface for the ODE/DAE solvers is TSSolve(), you should almost for sure be using that routine and not this routine. 3029 3030 The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may 3031 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 3032 3033 This may over-step the final time provided in TSSetDuration() depending on the time-step used. TSSolve() interpolates to exactly the 3034 time provided in TSSetDuration(). One can use TSInterpolate() to determine an interpolated solution within the final timestep. 3035 3036 .keywords: TS, timestep, solve 3037 3038 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSInterpolate() 3039 @*/ 3040 PetscErrorCode TSStep(TS ts) 3041 { 3042 DM dm; 3043 PetscErrorCode ierr; 3044 static PetscBool cite = PETSC_FALSE; 3045 3046 PetscFunctionBegin; 3047 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 3048 ierr = PetscCitationsRegister("@techreport{tspaper,\n" 3049 " title = {{PETSc/TS}: A Modern Scalable {DAE/ODE} Solver Library},\n" 3050 " author = {Shrirang Abhyankar and Jed Brown and Emil Constantinescu and Debojyoti Ghosh and Barry F. Smith},\n" 3051 " type = {Preprint},\n" 3052 " number = {ANL/MCS-P5061-0114},\n" 3053 " institution = {Argonne National Laboratory},\n" 3054 " year = {2014}\n}\n",&cite); 3055 3056 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 3057 ierr = TSSetUp(ts);CHKERRQ(ierr); 3058 3059 ts->reason = TS_CONVERGED_ITERATING; 3060 ts->ptime_prev = ts->ptime; 3061 ierr = DMSetOutputSequenceNumber(dm, ts->steps, ts->ptime);CHKERRQ(ierr); 3062 3063 if (!ts->ops->step) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSStep not implemented for type '%s'",((PetscObject)ts)->type_name); 3064 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 3065 ierr = (*ts->ops->step)(ts);CHKERRQ(ierr); 3066 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 3067 3068 ts->time_step_prev = ts->ptime - ts->ptime_prev; 3069 ierr = DMSetOutputSequenceNumber(dm, ts->steps, ts->ptime);CHKERRQ(ierr); 3070 3071 if (ts->reason < 0) { 3072 if (ts->errorifstepfailed) { 3073 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]); 3074 else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 3075 } 3076 } else if (!ts->reason) { 3077 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 3078 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 3079 } 3080 ts->total_steps++; 3081 PetscFunctionReturn(0); 3082 } 3083 3084 #undef __FUNCT__ 3085 #define __FUNCT__ "TSAdjointStep" 3086 /*@ 3087 TSAdjointStep - Steps one time step 3088 3089 Collective on TS 3090 3091 Input Parameter: 3092 . ts - the TS context obtained from TSCreate() 3093 3094 Level: intermediate 3095 3096 Notes: 3097 The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may 3098 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 3099 3100 This may over-step the final time provided in TSSetDuration() depending on the time-step used. TSSolve() interpolates to exactly the 3101 time provided in TSSetDuration(). One can use TSInterpolate() to determine an interpolated solution within the final timestep. 3102 3103 .keywords: TS, timestep, solve 3104 3105 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSInterpolate() 3106 @*/ 3107 PetscErrorCode TSAdjointStep(TS ts) 3108 { 3109 DM dm; 3110 PetscErrorCode ierr; 3111 3112 PetscFunctionBegin; 3113 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 3114 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 3115 ierr = TSAdjointSetUp(ts);CHKERRQ(ierr); 3116 3117 ts->reason = TS_CONVERGED_ITERATING; 3118 ts->ptime_prev = ts->ptime; 3119 ierr = DMSetOutputSequenceNumber(dm, ts->steps, ts->ptime);CHKERRQ(ierr); 3120 ierr = VecViewFromOptions(ts->vec_sol, ((PetscObject) ts)->prefix, "-ts_view_solution");CHKERRQ(ierr); 3121 3122 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 3123 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); 3124 ierr = (*ts->ops->adjointstep)(ts);CHKERRQ(ierr); 3125 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 3126 3127 ts->time_step_prev = ts->ptime - ts->ptime_prev; 3128 ierr = DMSetOutputSequenceNumber(dm, ts->steps, ts->ptime);CHKERRQ(ierr); 3129 3130 if (ts->reason < 0) { 3131 if (ts->errorifstepfailed) { 3132 if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) { 3133 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]); 3134 } else if (ts->reason == TS_DIVERGED_STEP_REJECTED) { 3135 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]); 3136 } else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 3137 } 3138 } else if (!ts->reason) { 3139 if (ts->steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS; 3140 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 3141 } 3142 ts->total_steps--; 3143 PetscFunctionReturn(0); 3144 } 3145 3146 #undef __FUNCT__ 3147 #define __FUNCT__ "TSEvaluateStep" 3148 /*@ 3149 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 3150 3151 Collective on TS 3152 3153 Input Arguments: 3154 + ts - time stepping context 3155 . order - desired order of accuracy 3156 - done - whether the step was evaluated at this order (pass NULL to generate an error if not available) 3157 3158 Output Arguments: 3159 . U - state at the end of the current step 3160 3161 Level: advanced 3162 3163 Notes: 3164 This function cannot be called until all stages have been evaluated. 3165 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. 3166 3167 .seealso: TSStep(), TSAdapt 3168 @*/ 3169 PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done) 3170 { 3171 PetscErrorCode ierr; 3172 3173 PetscFunctionBegin; 3174 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3175 PetscValidType(ts,1); 3176 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 3177 if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name); 3178 ierr = (*ts->ops->evaluatestep)(ts,order,U,done);CHKERRQ(ierr); 3179 PetscFunctionReturn(0); 3180 } 3181 3182 3183 #undef __FUNCT__ 3184 #define __FUNCT__ "TSSolve" 3185 /*@ 3186 TSSolve - Steps the requested number of timesteps. 3187 3188 Collective on TS 3189 3190 Input Parameter: 3191 + ts - the TS context obtained from TSCreate() 3192 - u - the solution vector (can be null if TSSetSolution() was used, otherwise must contain the initial conditions) 3193 3194 Level: beginner 3195 3196 Notes: 3197 The final time returned by this function may be different from the time of the internally 3198 held state accessible by TSGetSolution() and TSGetTime() because the method may have 3199 stepped over the final time. 3200 3201 .keywords: TS, timestep, solve 3202 3203 .seealso: TSCreate(), TSSetSolution(), TSStep() 3204 @*/ 3205 PetscErrorCode TSSolve(TS ts,Vec u) 3206 { 3207 Vec solution; 3208 PetscErrorCode ierr; 3209 3210 PetscFunctionBegin; 3211 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3212 if (u) PetscValidHeaderSpecific(u,VEC_CLASSID,2); 3213 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 */ 3214 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 3215 if (!ts->vec_sol || u == ts->vec_sol) { 3216 ierr = VecDuplicate(u,&solution);CHKERRQ(ierr); 3217 ierr = TSSetSolution(ts,solution);CHKERRQ(ierr); 3218 ierr = VecDestroy(&solution);CHKERRQ(ierr); /* grant ownership */ 3219 } 3220 ierr = VecCopy(u,ts->vec_sol);CHKERRQ(ierr); 3221 } else if (u) { 3222 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 3223 } 3224 ierr = TSSetUp(ts);CHKERRQ(ierr); /*compute adj coefficients if the reverse mode is on*/ 3225 /* reset time step and iteration counters */ 3226 ts->steps = 0; 3227 ts->ksp_its = 0; 3228 ts->snes_its = 0; 3229 ts->num_snes_failures = 0; 3230 ts->reject = 0; 3231 ts->reason = TS_CONVERGED_ITERATING; 3232 3233 ierr = TSViewFromOptions(ts,NULL,"-ts_view_pre");CHKERRQ(ierr); 3234 3235 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 3236 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 3237 ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr); 3238 ts->solvetime = ts->ptime; 3239 } else { 3240 /* steps the requested number of timesteps. */ 3241 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 3242 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 3243 while (!ts->reason) { 3244 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 3245 ierr = TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 3246 ierr = TSStep(ts);CHKERRQ(ierr); 3247 if (ts->event) { 3248 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 3249 if (ts->event->status != TSEVENT_PROCESSING) { 3250 ierr = TSPostStep(ts);CHKERRQ(ierr); 3251 } 3252 } else { 3253 ierr = TSPostStep(ts);CHKERRQ(ierr); 3254 } 3255 } 3256 if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) { 3257 ierr = TSInterpolate(ts,ts->max_time,u);CHKERRQ(ierr); 3258 ts->solvetime = ts->max_time; 3259 solution = u; 3260 } else { 3261 if (u) {ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr);} 3262 ts->solvetime = ts->ptime; 3263 solution = ts->vec_sol; 3264 } 3265 ierr = TSTrajectorySet(ts->trajectory,ts,ts->steps,ts->solvetime,solution);CHKERRQ(ierr); 3266 ierr = TSMonitor(ts,ts->steps,ts->solvetime,solution);CHKERRQ(ierr); 3267 ierr = VecViewFromOptions(solution, ((PetscObject) ts)->prefix, "-ts_view_solution");CHKERRQ(ierr); 3268 } 3269 3270 ierr = TSViewFromOptions(ts,NULL,"-ts_view");CHKERRQ(ierr); 3271 ierr = PetscObjectSAWsBlock((PetscObject)ts);CHKERRQ(ierr); 3272 PetscFunctionReturn(0); 3273 } 3274 3275 #undef __FUNCT__ 3276 #define __FUNCT__ "TSAdjointSolve" 3277 /*@ 3278 TSAdjointSolve - Solves the discrete ajoint problem for an ODE/DAE 3279 3280 Collective on TS 3281 3282 Input Parameter: 3283 . ts - the TS context obtained from TSCreate() 3284 3285 Level: intermediate 3286 3287 Notes: 3288 This must be called after a call to TSSolve() that solves the forward problem 3289 3290 By default this will integrate back to the initial time, one can use TSAdjointSetSteps() to step back to a later time 3291 3292 .keywords: TS, timestep, solve 3293 3294 .seealso: TSCreate(), TSSetSolution(), TSStep() 3295 @*/ 3296 PetscErrorCode TSAdjointSolve(TS ts) 3297 { 3298 PetscErrorCode ierr; 3299 3300 PetscFunctionBegin; 3301 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3302 ierr = TSAdjointSetUp(ts);CHKERRQ(ierr); 3303 /* reset time step and iteration counters */ 3304 ts->steps = 0; 3305 ts->ksp_its = 0; 3306 ts->snes_its = 0; 3307 ts->num_snes_failures = 0; 3308 ts->reject = 0; 3309 ts->reason = TS_CONVERGED_ITERATING; 3310 3311 if (!ts->adjoint_max_steps) ts->adjoint_max_steps = ts->total_steps; 3312 3313 if (ts->steps >= ts->adjoint_max_steps) ts->reason = TS_CONVERGED_ITS; 3314 while (!ts->reason) { 3315 ierr = TSTrajectoryGet(ts->trajectory,ts,ts->adjoint_max_steps-ts->steps,ts->ptime);CHKERRQ(ierr); 3316 ierr = TSMonitor(ts,ts->adjoint_max_steps-ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 3317 ierr = TSAdjointStep(ts);CHKERRQ(ierr); 3318 if (ts->event) { 3319 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 3320 if (ts->event->status != TSEVENT_PROCESSING) { 3321 ierr = TSPostStep(ts);CHKERRQ(ierr); 3322 } 3323 } else { 3324 ierr = TSPostStep(ts);CHKERRQ(ierr); 3325 } 3326 } 3327 ts->solvetime = ts->ptime; 3328 PetscFunctionReturn(0); 3329 } 3330 3331 #undef __FUNCT__ 3332 #define __FUNCT__ "TSMonitor" 3333 /*@ 3334 TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet() 3335 3336 Collective on TS 3337 3338 Input Parameters: 3339 + ts - time stepping context obtained from TSCreate() 3340 . step - step number that has just completed 3341 . ptime - model time of the state 3342 - u - state at the current model time 3343 3344 Notes: 3345 TSMonitor() is typically used within the time stepping implementations. 3346 Users might call this function when using the TSStep() interface instead of TSSolve(). 3347 3348 Level: advanced 3349 3350 .keywords: TS, timestep 3351 @*/ 3352 PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u) 3353 { 3354 PetscErrorCode ierr; 3355 PetscInt i,n = ts->numbermonitors; 3356 3357 PetscFunctionBegin; 3358 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3359 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 3360 ierr = VecLockPush(u);CHKERRQ(ierr); 3361 for (i=0; i<n; i++) { 3362 ierr = (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);CHKERRQ(ierr); 3363 } 3364 ierr = VecLockPop(u);CHKERRQ(ierr); 3365 PetscFunctionReturn(0); 3366 } 3367 3368 /* ------------------------------------------------------------------------*/ 3369 #undef __FUNCT__ 3370 #define __FUNCT__ "TSMonitorLGCtxCreate" 3371 /*@C 3372 TSMonitorLGCtxCreate - Creates a line graph context for use with 3373 TS to monitor the solution process graphically in various ways 3374 3375 Collective on TS 3376 3377 Input Parameters: 3378 + host - the X display to open, or null for the local machine 3379 . label - the title to put in the title bar 3380 . x, y - the screen coordinates of the upper left coordinate of the window 3381 . m, n - the screen width and height in pixels 3382 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time 3383 3384 Output Parameter: 3385 . ctx - the context 3386 3387 Options Database Key: 3388 + -ts_monitor_lg_timestep - automatically sets line graph monitor 3389 . -ts_monitor_lg_solution - 3390 . -ts_monitor_lg_error - 3391 . -ts_monitor_lg_ksp_iterations - 3392 . -ts_monitor_lg_snes_iterations - 3393 - -lg_indicate_data_points <true,false> - indicate the data points (at each time step) on the plot; default is true 3394 3395 Notes: 3396 Use TSMonitorLGCtxDestroy() to destroy. 3397 3398 Level: intermediate 3399 3400 .keywords: TS, monitor, line graph, residual, seealso 3401 3402 .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 3403 3404 @*/ 3405 PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx) 3406 { 3407 PetscDraw win; 3408 PetscErrorCode ierr; 3409 3410 PetscFunctionBegin; 3411 ierr = PetscNew(ctx);CHKERRQ(ierr); 3412 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 3413 ierr = PetscDrawSetFromOptions(win);CHKERRQ(ierr); 3414 ierr = PetscDrawLGCreate(win,1,&(*ctx)->lg);CHKERRQ(ierr); 3415 ierr = PetscLogObjectParent((PetscObject)(*ctx)->lg,(PetscObject)win);CHKERRQ(ierr); 3416 ierr = PetscDrawLGIndicateDataPoints((*ctx)->lg,PETSC_TRUE);CHKERRQ(ierr); 3417 ierr = PetscDrawLGSetFromOptions((*ctx)->lg);CHKERRQ(ierr); 3418 (*ctx)->howoften = howoften; 3419 PetscFunctionReturn(0); 3420 } 3421 3422 #undef __FUNCT__ 3423 #define __FUNCT__ "TSMonitorLGTimeStep" 3424 PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx) 3425 { 3426 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 3427 PetscReal x = ptime,y; 3428 PetscErrorCode ierr; 3429 3430 PetscFunctionBegin; 3431 if (!step) { 3432 PetscDrawAxis axis; 3433 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 3434 ierr = PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time","Time step");CHKERRQ(ierr); 3435 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 3436 ierr = PetscDrawLGIndicateDataPoints(ctx->lg,PETSC_TRUE);CHKERRQ(ierr); 3437 } 3438 ierr = TSGetTimeStep(ts,&y);CHKERRQ(ierr); 3439 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 3440 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 3441 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 3442 } 3443 PetscFunctionReturn(0); 3444 } 3445 3446 #undef __FUNCT__ 3447 #define __FUNCT__ "TSMonitorLGCtxDestroy" 3448 /*@C 3449 TSMonitorLGCtxDestroy - Destroys a line graph context that was created 3450 with TSMonitorLGCtxCreate(). 3451 3452 Collective on TSMonitorLGCtx 3453 3454 Input Parameter: 3455 . ctx - the monitor context 3456 3457 Level: intermediate 3458 3459 .keywords: TS, monitor, line graph, destroy 3460 3461 .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep(); 3462 @*/ 3463 PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx) 3464 { 3465 PetscDraw draw; 3466 PetscErrorCode ierr; 3467 3468 PetscFunctionBegin; 3469 ierr = PetscDrawLGGetDraw((*ctx)->lg,&draw);CHKERRQ(ierr); 3470 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 3471 ierr = PetscDrawLGDestroy(&(*ctx)->lg);CHKERRQ(ierr); 3472 ierr = PetscFree(*ctx);CHKERRQ(ierr); 3473 PetscFunctionReturn(0); 3474 } 3475 3476 #undef __FUNCT__ 3477 #define __FUNCT__ "TSGetTime" 3478 /*@ 3479 TSGetTime - Gets the time of the most recently completed step. 3480 3481 Not Collective 3482 3483 Input Parameter: 3484 . ts - the TS context obtained from TSCreate() 3485 3486 Output Parameter: 3487 . t - the current time 3488 3489 Level: beginner 3490 3491 Note: 3492 When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(), 3493 TSSetPreStage(), TSSetPostStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated. 3494 3495 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 3496 3497 .keywords: TS, get, time 3498 @*/ 3499 PetscErrorCode TSGetTime(TS ts,PetscReal *t) 3500 { 3501 PetscFunctionBegin; 3502 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3503 PetscValidRealPointer(t,2); 3504 *t = ts->ptime; 3505 PetscFunctionReturn(0); 3506 } 3507 3508 #undef __FUNCT__ 3509 #define __FUNCT__ "TSGetPrevTime" 3510 /*@ 3511 TSGetPrevTime - Gets the starting time of the previously completed step. 3512 3513 Not Collective 3514 3515 Input Parameter: 3516 . ts - the TS context obtained from TSCreate() 3517 3518 Output Parameter: 3519 . t - the previous time 3520 3521 Level: beginner 3522 3523 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 3524 3525 .keywords: TS, get, time 3526 @*/ 3527 PetscErrorCode TSGetPrevTime(TS ts,PetscReal *t) 3528 { 3529 PetscFunctionBegin; 3530 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3531 PetscValidRealPointer(t,2); 3532 *t = ts->ptime_prev; 3533 PetscFunctionReturn(0); 3534 } 3535 3536 #undef __FUNCT__ 3537 #define __FUNCT__ "TSSetTime" 3538 /*@ 3539 TSSetTime - Allows one to reset the time. 3540 3541 Logically Collective on TS 3542 3543 Input Parameters: 3544 + ts - the TS context obtained from TSCreate() 3545 - time - the time 3546 3547 Level: intermediate 3548 3549 .seealso: TSGetTime(), TSSetDuration() 3550 3551 .keywords: TS, set, time 3552 @*/ 3553 PetscErrorCode TSSetTime(TS ts, PetscReal t) 3554 { 3555 PetscFunctionBegin; 3556 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3557 PetscValidLogicalCollectiveReal(ts,t,2); 3558 ts->ptime = t; 3559 PetscFunctionReturn(0); 3560 } 3561 3562 #undef __FUNCT__ 3563 #define __FUNCT__ "TSSetOptionsPrefix" 3564 /*@C 3565 TSSetOptionsPrefix - Sets the prefix used for searching for all 3566 TS options in the database. 3567 3568 Logically Collective on TS 3569 3570 Input Parameter: 3571 + ts - The TS context 3572 - prefix - The prefix to prepend to all option names 3573 3574 Notes: 3575 A hyphen (-) must NOT be given at the beginning of the prefix name. 3576 The first character of all runtime options is AUTOMATICALLY the 3577 hyphen. 3578 3579 Level: advanced 3580 3581 .keywords: TS, set, options, prefix, database 3582 3583 .seealso: TSSetFromOptions() 3584 3585 @*/ 3586 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 3587 { 3588 PetscErrorCode ierr; 3589 SNES snes; 3590 3591 PetscFunctionBegin; 3592 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3593 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3594 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3595 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 3596 PetscFunctionReturn(0); 3597 } 3598 3599 3600 #undef __FUNCT__ 3601 #define __FUNCT__ "TSAppendOptionsPrefix" 3602 /*@C 3603 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 3604 TS options in the database. 3605 3606 Logically Collective on TS 3607 3608 Input Parameter: 3609 + ts - The TS context 3610 - prefix - The prefix to prepend to all option names 3611 3612 Notes: 3613 A hyphen (-) must NOT be given at the beginning of the prefix name. 3614 The first character of all runtime options is AUTOMATICALLY the 3615 hyphen. 3616 3617 Level: advanced 3618 3619 .keywords: TS, append, options, prefix, database 3620 3621 .seealso: TSGetOptionsPrefix() 3622 3623 @*/ 3624 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 3625 { 3626 PetscErrorCode ierr; 3627 SNES snes; 3628 3629 PetscFunctionBegin; 3630 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3631 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3632 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3633 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 3634 PetscFunctionReturn(0); 3635 } 3636 3637 #undef __FUNCT__ 3638 #define __FUNCT__ "TSGetOptionsPrefix" 3639 /*@C 3640 TSGetOptionsPrefix - Sets the prefix used for searching for all 3641 TS options in the database. 3642 3643 Not Collective 3644 3645 Input Parameter: 3646 . ts - The TS context 3647 3648 Output Parameter: 3649 . prefix - A pointer to the prefix string used 3650 3651 Notes: On the fortran side, the user should pass in a string 'prifix' of 3652 sufficient length to hold the prefix. 3653 3654 Level: intermediate 3655 3656 .keywords: TS, get, options, prefix, database 3657 3658 .seealso: TSAppendOptionsPrefix() 3659 @*/ 3660 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 3661 { 3662 PetscErrorCode ierr; 3663 3664 PetscFunctionBegin; 3665 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3666 PetscValidPointer(prefix,2); 3667 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 3668 PetscFunctionReturn(0); 3669 } 3670 3671 #undef __FUNCT__ 3672 #define __FUNCT__ "TSGetRHSJacobian" 3673 /*@C 3674 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 3675 3676 Not Collective, but parallel objects are returned if TS is parallel 3677 3678 Input Parameter: 3679 . ts - The TS context obtained from TSCreate() 3680 3681 Output Parameters: 3682 + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or NULL) 3683 . Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat (or NULL) 3684 . func - Function to compute the Jacobian of the RHS (or NULL) 3685 - ctx - User-defined context for Jacobian evaluation routine (or NULL) 3686 3687 Notes: You can pass in NULL for any return argument you do not need. 3688 3689 Level: intermediate 3690 3691 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3692 3693 .keywords: TS, timestep, get, matrix, Jacobian 3694 @*/ 3695 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx) 3696 { 3697 PetscErrorCode ierr; 3698 SNES snes; 3699 DM dm; 3700 3701 PetscFunctionBegin; 3702 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3703 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3704 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3705 ierr = DMTSGetRHSJacobian(dm,func,ctx);CHKERRQ(ierr); 3706 PetscFunctionReturn(0); 3707 } 3708 3709 #undef __FUNCT__ 3710 #define __FUNCT__ "TSGetIJacobian" 3711 /*@C 3712 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 3713 3714 Not Collective, but parallel objects are returned if TS is parallel 3715 3716 Input Parameter: 3717 . ts - The TS context obtained from TSCreate() 3718 3719 Output Parameters: 3720 + Amat - The (approximate) Jacobian of F(t,U,U_t) 3721 . Pmat - The matrix from which the preconditioner is constructed, often the same as Amat 3722 . f - The function to compute the matrices 3723 - ctx - User-defined context for Jacobian evaluation routine 3724 3725 Notes: You can pass in NULL for any return argument you do not need. 3726 3727 Level: advanced 3728 3729 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3730 3731 .keywords: TS, timestep, get, matrix, Jacobian 3732 @*/ 3733 PetscErrorCode TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx) 3734 { 3735 PetscErrorCode ierr; 3736 SNES snes; 3737 DM dm; 3738 3739 PetscFunctionBegin; 3740 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3741 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 3742 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3743 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3744 ierr = DMTSGetIJacobian(dm,f,ctx);CHKERRQ(ierr); 3745 PetscFunctionReturn(0); 3746 } 3747 3748 3749 #undef __FUNCT__ 3750 #define __FUNCT__ "TSMonitorDrawSolution" 3751 /*@C 3752 TSMonitorDrawSolution - Monitors progress of the TS solvers by calling 3753 VecView() for the solution at each timestep 3754 3755 Collective on TS 3756 3757 Input Parameters: 3758 + ts - the TS context 3759 . step - current time-step 3760 . ptime - current time 3761 - dummy - either a viewer or NULL 3762 3763 Options Database: 3764 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3765 3766 Notes: the initial solution and current solution are not displayed with a common axis scaling so generally the option -ts_monitor_draw_solution_initial 3767 will look bad 3768 3769 Level: intermediate 3770 3771 .keywords: TS, vector, monitor, view 3772 3773 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3774 @*/ 3775 PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3776 { 3777 PetscErrorCode ierr; 3778 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3779 PetscDraw draw; 3780 3781 PetscFunctionBegin; 3782 if (!step && ictx->showinitial) { 3783 if (!ictx->initialsolution) { 3784 ierr = VecDuplicate(u,&ictx->initialsolution);CHKERRQ(ierr); 3785 } 3786 ierr = VecCopy(u,ictx->initialsolution);CHKERRQ(ierr); 3787 } 3788 if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 3789 3790 if (ictx->showinitial) { 3791 PetscReal pause; 3792 ierr = PetscViewerDrawGetPause(ictx->viewer,&pause);CHKERRQ(ierr); 3793 ierr = PetscViewerDrawSetPause(ictx->viewer,0.0);CHKERRQ(ierr); 3794 ierr = VecView(ictx->initialsolution,ictx->viewer);CHKERRQ(ierr); 3795 ierr = PetscViewerDrawSetPause(ictx->viewer,pause);CHKERRQ(ierr); 3796 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);CHKERRQ(ierr); 3797 } 3798 ierr = VecView(u,ictx->viewer);CHKERRQ(ierr); 3799 if (ictx->showtimestepandtime) { 3800 PetscReal xl,yl,xr,yr,tw,w,h; 3801 char time[32]; 3802 size_t len; 3803 3804 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3805 ierr = PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);CHKERRQ(ierr); 3806 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3807 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3808 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3809 w = xl + .5*(xr - xl) - .5*len*tw; 3810 h = yl + .95*(yr - yl); 3811 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3812 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3813 } 3814 3815 if (ictx->showinitial) { 3816 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);CHKERRQ(ierr); 3817 } 3818 PetscFunctionReturn(0); 3819 } 3820 3821 #undef __FUNCT__ 3822 #define __FUNCT__ "TSMonitorDrawSolutionPhase" 3823 /*@C 3824 TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram 3825 3826 Collective on TS 3827 3828 Input Parameters: 3829 + ts - the TS context 3830 . step - current time-step 3831 . ptime - current time 3832 - dummy - either a viewer or NULL 3833 3834 Level: intermediate 3835 3836 .keywords: TS, vector, monitor, view 3837 3838 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3839 @*/ 3840 PetscErrorCode TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3841 { 3842 PetscErrorCode ierr; 3843 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3844 PetscDraw draw; 3845 MPI_Comm comm; 3846 PetscInt n; 3847 PetscMPIInt size; 3848 PetscReal xl,yl,xr,yr,tw,w,h; 3849 char time[32]; 3850 size_t len; 3851 const PetscScalar *U; 3852 3853 PetscFunctionBegin; 3854 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 3855 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3856 if (size != 1) SETERRQ(comm,PETSC_ERR_SUP,"Only allowed for sequential runs"); 3857 ierr = VecGetSize(u,&n);CHKERRQ(ierr); 3858 if (n != 2) SETERRQ(comm,PETSC_ERR_SUP,"Only for ODEs with two unknowns"); 3859 3860 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3861 3862 ierr = VecGetArrayRead(u,&U);CHKERRQ(ierr); 3863 ierr = PetscDrawAxisGetLimits(ictx->axis,&xl,&xr,&yl,&yr);CHKERRQ(ierr); 3864 if ((PetscRealPart(U[0]) < xl) || (PetscRealPart(U[1]) < yl) || (PetscRealPart(U[0]) > xr) || (PetscRealPart(U[1]) > yr)) { 3865 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3866 PetscFunctionReturn(0); 3867 } 3868 if (!step) ictx->color++; 3869 ierr = PetscDrawPoint(draw,PetscRealPart(U[0]),PetscRealPart(U[1]),ictx->color);CHKERRQ(ierr); 3870 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3871 3872 if (ictx->showtimestepandtime) { 3873 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3874 ierr = PetscSNPrintf(time,32,"Timestep %d Time %g",(int)step,(double)ptime);CHKERRQ(ierr); 3875 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3876 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3877 w = xl + .5*(xr - xl) - .5*len*tw; 3878 h = yl + .95*(yr - yl); 3879 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3880 } 3881 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3882 PetscFunctionReturn(0); 3883 } 3884 3885 3886 #undef __FUNCT__ 3887 #define __FUNCT__ "TSMonitorDrawCtxDestroy" 3888 /*@C 3889 TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution() 3890 3891 Collective on TS 3892 3893 Input Parameters: 3894 . ctx - the monitor context 3895 3896 Level: intermediate 3897 3898 .keywords: TS, vector, monitor, view 3899 3900 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError() 3901 @*/ 3902 PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx) 3903 { 3904 PetscErrorCode ierr; 3905 3906 PetscFunctionBegin; 3907 ierr = PetscDrawAxisDestroy(&(*ictx)->axis);CHKERRQ(ierr); 3908 ierr = PetscViewerDestroy(&(*ictx)->viewer);CHKERRQ(ierr); 3909 ierr = VecDestroy(&(*ictx)->initialsolution);CHKERRQ(ierr); 3910 ierr = PetscFree(*ictx);CHKERRQ(ierr); 3911 PetscFunctionReturn(0); 3912 } 3913 3914 #undef __FUNCT__ 3915 #define __FUNCT__ "TSMonitorDrawCtxCreate" 3916 /*@C 3917 TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx 3918 3919 Collective on TS 3920 3921 Input Parameter: 3922 . ts - time-step context 3923 3924 Output Patameter: 3925 . ctx - the monitor context 3926 3927 Options Database: 3928 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3929 3930 Level: intermediate 3931 3932 .keywords: TS, vector, monitor, view 3933 3934 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx() 3935 @*/ 3936 PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx) 3937 { 3938 PetscErrorCode ierr; 3939 3940 PetscFunctionBegin; 3941 ierr = PetscNew(ctx);CHKERRQ(ierr); 3942 ierr = PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);CHKERRQ(ierr); 3943 ierr = PetscViewerSetFromOptions((*ctx)->viewer);CHKERRQ(ierr); 3944 3945 (*ctx)->howoften = howoften; 3946 (*ctx)->showinitial = PETSC_FALSE; 3947 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);CHKERRQ(ierr); 3948 3949 (*ctx)->showtimestepandtime = PETSC_FALSE; 3950 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);CHKERRQ(ierr); 3951 (*ctx)->color = PETSC_DRAW_WHITE; 3952 PetscFunctionReturn(0); 3953 } 3954 3955 #undef __FUNCT__ 3956 #define __FUNCT__ "TSMonitorDrawError" 3957 /*@C 3958 TSMonitorDrawError - Monitors progress of the TS solvers by calling 3959 VecView() for the error at each timestep 3960 3961 Collective on TS 3962 3963 Input Parameters: 3964 + ts - the TS context 3965 . step - current time-step 3966 . ptime - current time 3967 - dummy - either a viewer or NULL 3968 3969 Level: intermediate 3970 3971 .keywords: TS, vector, monitor, view 3972 3973 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3974 @*/ 3975 PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3976 { 3977 PetscErrorCode ierr; 3978 TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy; 3979 PetscViewer viewer = ctx->viewer; 3980 Vec work; 3981 3982 PetscFunctionBegin; 3983 if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 3984 ierr = VecDuplicate(u,&work);CHKERRQ(ierr); 3985 ierr = TSComputeSolutionFunction(ts,ptime,work);CHKERRQ(ierr); 3986 ierr = VecAXPY(work,-1.0,u);CHKERRQ(ierr); 3987 ierr = VecView(work,viewer);CHKERRQ(ierr); 3988 ierr = VecDestroy(&work);CHKERRQ(ierr); 3989 PetscFunctionReturn(0); 3990 } 3991 3992 #include <petsc-private/dmimpl.h> 3993 #undef __FUNCT__ 3994 #define __FUNCT__ "TSSetDM" 3995 /*@ 3996 TSSetDM - Sets the DM that may be used by some preconditioners 3997 3998 Logically Collective on TS and DM 3999 4000 Input Parameters: 4001 + ts - the preconditioner context 4002 - dm - the dm 4003 4004 Level: intermediate 4005 4006 4007 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 4008 @*/ 4009 PetscErrorCode TSSetDM(TS ts,DM dm) 4010 { 4011 PetscErrorCode ierr; 4012 SNES snes; 4013 DMTS tsdm; 4014 4015 PetscFunctionBegin; 4016 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4017 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 4018 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 4019 if (ts->dm->dmts && !dm->dmts) { 4020 ierr = DMCopyDMTS(ts->dm,dm);CHKERRQ(ierr); 4021 ierr = DMGetDMTS(ts->dm,&tsdm);CHKERRQ(ierr); 4022 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 4023 tsdm->originaldm = dm; 4024 } 4025 } 4026 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 4027 } 4028 ts->dm = dm; 4029 4030 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4031 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 4032 PetscFunctionReturn(0); 4033 } 4034 4035 #undef __FUNCT__ 4036 #define __FUNCT__ "TSGetDM" 4037 /*@ 4038 TSGetDM - Gets the DM that may be used by some preconditioners 4039 4040 Not Collective 4041 4042 Input Parameter: 4043 . ts - the preconditioner context 4044 4045 Output Parameter: 4046 . dm - the dm 4047 4048 Level: intermediate 4049 4050 4051 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 4052 @*/ 4053 PetscErrorCode TSGetDM(TS ts,DM *dm) 4054 { 4055 PetscErrorCode ierr; 4056 4057 PetscFunctionBegin; 4058 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4059 if (!ts->dm) { 4060 ierr = DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);CHKERRQ(ierr); 4061 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 4062 } 4063 *dm = ts->dm; 4064 PetscFunctionReturn(0); 4065 } 4066 4067 #undef __FUNCT__ 4068 #define __FUNCT__ "SNESTSFormFunction" 4069 /*@ 4070 SNESTSFormFunction - Function to evaluate nonlinear residual 4071 4072 Logically Collective on SNES 4073 4074 Input Parameter: 4075 + snes - nonlinear solver 4076 . U - the current state at which to evaluate the residual 4077 - ctx - user context, must be a TS 4078 4079 Output Parameter: 4080 . F - the nonlinear residual 4081 4082 Notes: 4083 This function is not normally called by users and is automatically registered with the SNES used by TS. 4084 It is most frequently passed to MatFDColoringSetFunction(). 4085 4086 Level: advanced 4087 4088 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 4089 @*/ 4090 PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx) 4091 { 4092 TS ts = (TS)ctx; 4093 PetscErrorCode ierr; 4094 4095 PetscFunctionBegin; 4096 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4097 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 4098 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 4099 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 4100 ierr = (ts->ops->snesfunction)(snes,U,F,ts);CHKERRQ(ierr); 4101 PetscFunctionReturn(0); 4102 } 4103 4104 #undef __FUNCT__ 4105 #define __FUNCT__ "SNESTSFormJacobian" 4106 /*@ 4107 SNESTSFormJacobian - Function to evaluate the Jacobian 4108 4109 Collective on SNES 4110 4111 Input Parameter: 4112 + snes - nonlinear solver 4113 . U - the current state at which to evaluate the residual 4114 - ctx - user context, must be a TS 4115 4116 Output Parameter: 4117 + A - the Jacobian 4118 . B - the preconditioning matrix (may be the same as A) 4119 - flag - indicates any structure change in the matrix 4120 4121 Notes: 4122 This function is not normally called by users and is automatically registered with the SNES used by TS. 4123 4124 Level: developer 4125 4126 .seealso: SNESSetJacobian() 4127 @*/ 4128 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat A,Mat B,void *ctx) 4129 { 4130 TS ts = (TS)ctx; 4131 PetscErrorCode ierr; 4132 4133 PetscFunctionBegin; 4134 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4135 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 4136 PetscValidPointer(A,3); 4137 PetscValidHeaderSpecific(A,MAT_CLASSID,3); 4138 PetscValidPointer(B,4); 4139 PetscValidHeaderSpecific(B,MAT_CLASSID,4); 4140 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 4141 ierr = (ts->ops->snesjacobian)(snes,U,A,B,ts);CHKERRQ(ierr); 4142 PetscFunctionReturn(0); 4143 } 4144 4145 #undef __FUNCT__ 4146 #define __FUNCT__ "TSComputeRHSFunctionLinear" 4147 /*@C 4148 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 4149 4150 Collective on TS 4151 4152 Input Arguments: 4153 + ts - time stepping context 4154 . t - time at which to evaluate 4155 . U - state at which to evaluate 4156 - ctx - context 4157 4158 Output Arguments: 4159 . F - right hand side 4160 4161 Level: intermediate 4162 4163 Notes: 4164 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 4165 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 4166 4167 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 4168 @*/ 4169 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 4170 { 4171 PetscErrorCode ierr; 4172 Mat Arhs,Brhs; 4173 4174 PetscFunctionBegin; 4175 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 4176 ierr = TSComputeRHSJacobian(ts,t,U,Arhs,Brhs);CHKERRQ(ierr); 4177 ierr = MatMult(Arhs,U,F);CHKERRQ(ierr); 4178 PetscFunctionReturn(0); 4179 } 4180 4181 #undef __FUNCT__ 4182 #define __FUNCT__ "TSComputeRHSJacobianConstant" 4183 /*@C 4184 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 4185 4186 Collective on TS 4187 4188 Input Arguments: 4189 + ts - time stepping context 4190 . t - time at which to evaluate 4191 . U - state at which to evaluate 4192 - ctx - context 4193 4194 Output Arguments: 4195 + A - pointer to operator 4196 . B - pointer to preconditioning matrix 4197 - flg - matrix structure flag 4198 4199 Level: intermediate 4200 4201 Notes: 4202 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 4203 4204 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 4205 @*/ 4206 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx) 4207 { 4208 PetscFunctionBegin; 4209 PetscFunctionReturn(0); 4210 } 4211 4212 #undef __FUNCT__ 4213 #define __FUNCT__ "TSComputeIFunctionLinear" 4214 /*@C 4215 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 4216 4217 Collective on TS 4218 4219 Input Arguments: 4220 + ts - time stepping context 4221 . t - time at which to evaluate 4222 . U - state at which to evaluate 4223 . Udot - time derivative of state vector 4224 - ctx - context 4225 4226 Output Arguments: 4227 . F - left hand side 4228 4229 Level: intermediate 4230 4231 Notes: 4232 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 4233 user is required to write their own TSComputeIFunction. 4234 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 4235 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 4236 4237 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 4238 @*/ 4239 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 4240 { 4241 PetscErrorCode ierr; 4242 Mat A,B; 4243 4244 PetscFunctionBegin; 4245 ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr); 4246 ierr = TSComputeIJacobian(ts,t,U,Udot,1.0,A,B,PETSC_TRUE);CHKERRQ(ierr); 4247 ierr = MatMult(A,Udot,F);CHKERRQ(ierr); 4248 PetscFunctionReturn(0); 4249 } 4250 4251 #undef __FUNCT__ 4252 #define __FUNCT__ "TSComputeIJacobianConstant" 4253 /*@C 4254 TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE 4255 4256 Collective on TS 4257 4258 Input Arguments: 4259 + ts - time stepping context 4260 . t - time at which to evaluate 4261 . U - state at which to evaluate 4262 . Udot - time derivative of state vector 4263 . shift - shift to apply 4264 - ctx - context 4265 4266 Output Arguments: 4267 + A - pointer to operator 4268 . B - pointer to preconditioning matrix 4269 - flg - matrix structure flag 4270 4271 Level: advanced 4272 4273 Notes: 4274 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 4275 4276 It is only appropriate for problems of the form 4277 4278 $ M Udot = F(U,t) 4279 4280 where M is constant and F is non-stiff. The user must pass M to TSSetIJacobian(). The current implementation only 4281 works with IMEX time integration methods such as TSROSW and TSARKIMEX, since there is no support for de-constructing 4282 an implicit operator of the form 4283 4284 $ shift*M + J 4285 4286 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 4287 a copy of M or reassemble it when requested. 4288 4289 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 4290 @*/ 4291 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat A,Mat B,void *ctx) 4292 { 4293 PetscErrorCode ierr; 4294 4295 PetscFunctionBegin; 4296 ierr = MatScale(A, shift / ts->ijacobian.shift);CHKERRQ(ierr); 4297 ts->ijacobian.shift = shift; 4298 PetscFunctionReturn(0); 4299 } 4300 4301 #undef __FUNCT__ 4302 #define __FUNCT__ "TSGetEquationType" 4303 /*@ 4304 TSGetEquationType - Gets the type of the equation that TS is solving. 4305 4306 Not Collective 4307 4308 Input Parameter: 4309 . ts - the TS context 4310 4311 Output Parameter: 4312 . equation_type - see TSEquationType 4313 4314 Level: beginner 4315 4316 .keywords: TS, equation type 4317 4318 .seealso: TSSetEquationType(), TSEquationType 4319 @*/ 4320 PetscErrorCode TSGetEquationType(TS ts,TSEquationType *equation_type) 4321 { 4322 PetscFunctionBegin; 4323 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4324 PetscValidPointer(equation_type,2); 4325 *equation_type = ts->equation_type; 4326 PetscFunctionReturn(0); 4327 } 4328 4329 #undef __FUNCT__ 4330 #define __FUNCT__ "TSSetEquationType" 4331 /*@ 4332 TSSetEquationType - Sets the type of the equation that TS is solving. 4333 4334 Not Collective 4335 4336 Input Parameter: 4337 + ts - the TS context 4338 . equation_type - see TSEquationType 4339 4340 Level: advanced 4341 4342 .keywords: TS, equation type 4343 4344 .seealso: TSGetEquationType(), TSEquationType 4345 @*/ 4346 PetscErrorCode TSSetEquationType(TS ts,TSEquationType equation_type) 4347 { 4348 PetscFunctionBegin; 4349 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4350 ts->equation_type = equation_type; 4351 PetscFunctionReturn(0); 4352 } 4353 4354 #undef __FUNCT__ 4355 #define __FUNCT__ "TSGetConvergedReason" 4356 /*@ 4357 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 4358 4359 Not Collective 4360 4361 Input Parameter: 4362 . ts - the TS context 4363 4364 Output Parameter: 4365 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 4366 manual pages for the individual convergence tests for complete lists 4367 4368 Level: beginner 4369 4370 Notes: 4371 Can only be called after the call to TSSolve() is complete. 4372 4373 .keywords: TS, nonlinear, set, convergence, test 4374 4375 .seealso: TSSetConvergenceTest(), TSConvergedReason 4376 @*/ 4377 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 4378 { 4379 PetscFunctionBegin; 4380 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4381 PetscValidPointer(reason,2); 4382 *reason = ts->reason; 4383 PetscFunctionReturn(0); 4384 } 4385 4386 #undef __FUNCT__ 4387 #define __FUNCT__ "TSSetConvergedReason" 4388 /*@ 4389 TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve. 4390 4391 Not Collective 4392 4393 Input Parameter: 4394 + ts - the TS context 4395 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 4396 manual pages for the individual convergence tests for complete lists 4397 4398 Level: advanced 4399 4400 Notes: 4401 Can only be called during TSSolve() is active. 4402 4403 .keywords: TS, nonlinear, set, convergence, test 4404 4405 .seealso: TSConvergedReason 4406 @*/ 4407 PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason) 4408 { 4409 PetscFunctionBegin; 4410 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4411 ts->reason = reason; 4412 PetscFunctionReturn(0); 4413 } 4414 4415 #undef __FUNCT__ 4416 #define __FUNCT__ "TSGetSolveTime" 4417 /*@ 4418 TSGetSolveTime - Gets the time after a call to TSSolve() 4419 4420 Not Collective 4421 4422 Input Parameter: 4423 . ts - the TS context 4424 4425 Output Parameter: 4426 . ftime - the final time. This time should correspond to the final time set with TSSetDuration() 4427 4428 Level: beginner 4429 4430 Notes: 4431 Can only be called after the call to TSSolve() is complete. 4432 4433 .keywords: TS, nonlinear, set, convergence, test 4434 4435 .seealso: TSSetConvergenceTest(), TSConvergedReason 4436 @*/ 4437 PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime) 4438 { 4439 PetscFunctionBegin; 4440 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4441 PetscValidPointer(ftime,2); 4442 *ftime = ts->solvetime; 4443 PetscFunctionReturn(0); 4444 } 4445 4446 #undef __FUNCT__ 4447 #define __FUNCT__ "TSGetTotalSteps" 4448 /*@ 4449 TSGetTotalSteps - Gets the total number of steps done since the last call to TSSetUp() or TSCreate() 4450 4451 Not Collective 4452 4453 Input Parameter: 4454 . ts - the TS context 4455 4456 Output Parameter: 4457 . steps - the number of steps 4458 4459 Level: beginner 4460 4461 Notes: 4462 Includes the number of steps for all calls to TSSolve() since TSSetUp() was called 4463 4464 .keywords: TS, nonlinear, set, convergence, test 4465 4466 .seealso: TSSetConvergenceTest(), TSConvergedReason 4467 @*/ 4468 PetscErrorCode TSGetTotalSteps(TS ts,PetscInt *steps) 4469 { 4470 PetscFunctionBegin; 4471 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4472 PetscValidPointer(steps,2); 4473 *steps = ts->total_steps; 4474 PetscFunctionReturn(0); 4475 } 4476 4477 #undef __FUNCT__ 4478 #define __FUNCT__ "TSGetSNESIterations" 4479 /*@ 4480 TSGetSNESIterations - Gets the total number of nonlinear iterations 4481 used by the time integrator. 4482 4483 Not Collective 4484 4485 Input Parameter: 4486 . ts - TS context 4487 4488 Output Parameter: 4489 . nits - number of nonlinear iterations 4490 4491 Notes: 4492 This counter is reset to zero for each successive call to TSSolve(). 4493 4494 Level: intermediate 4495 4496 .keywords: TS, get, number, nonlinear, iterations 4497 4498 .seealso: TSGetKSPIterations() 4499 @*/ 4500 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 4501 { 4502 PetscFunctionBegin; 4503 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4504 PetscValidIntPointer(nits,2); 4505 *nits = ts->snes_its; 4506 PetscFunctionReturn(0); 4507 } 4508 4509 #undef __FUNCT__ 4510 #define __FUNCT__ "TSGetKSPIterations" 4511 /*@ 4512 TSGetKSPIterations - Gets the total number of linear iterations 4513 used by the time integrator. 4514 4515 Not Collective 4516 4517 Input Parameter: 4518 . ts - TS context 4519 4520 Output Parameter: 4521 . lits - number of linear iterations 4522 4523 Notes: 4524 This counter is reset to zero for each successive call to TSSolve(). 4525 4526 Level: intermediate 4527 4528 .keywords: TS, get, number, linear, iterations 4529 4530 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 4531 @*/ 4532 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 4533 { 4534 PetscFunctionBegin; 4535 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4536 PetscValidIntPointer(lits,2); 4537 *lits = ts->ksp_its; 4538 PetscFunctionReturn(0); 4539 } 4540 4541 #undef __FUNCT__ 4542 #define __FUNCT__ "TSGetStepRejections" 4543 /*@ 4544 TSGetStepRejections - Gets the total number of rejected steps. 4545 4546 Not Collective 4547 4548 Input Parameter: 4549 . ts - TS context 4550 4551 Output Parameter: 4552 . rejects - number of steps rejected 4553 4554 Notes: 4555 This counter is reset to zero for each successive call to TSSolve(). 4556 4557 Level: intermediate 4558 4559 .keywords: TS, get, number 4560 4561 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 4562 @*/ 4563 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 4564 { 4565 PetscFunctionBegin; 4566 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4567 PetscValidIntPointer(rejects,2); 4568 *rejects = ts->reject; 4569 PetscFunctionReturn(0); 4570 } 4571 4572 #undef __FUNCT__ 4573 #define __FUNCT__ "TSGetSNESFailures" 4574 /*@ 4575 TSGetSNESFailures - Gets the total number of failed SNES solves 4576 4577 Not Collective 4578 4579 Input Parameter: 4580 . ts - TS context 4581 4582 Output Parameter: 4583 . fails - number of failed nonlinear solves 4584 4585 Notes: 4586 This counter is reset to zero for each successive call to TSSolve(). 4587 4588 Level: intermediate 4589 4590 .keywords: TS, get, number 4591 4592 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 4593 @*/ 4594 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 4595 { 4596 PetscFunctionBegin; 4597 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4598 PetscValidIntPointer(fails,2); 4599 *fails = ts->num_snes_failures; 4600 PetscFunctionReturn(0); 4601 } 4602 4603 #undef __FUNCT__ 4604 #define __FUNCT__ "TSSetMaxStepRejections" 4605 /*@ 4606 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 4607 4608 Not Collective 4609 4610 Input Parameter: 4611 + ts - TS context 4612 - rejects - maximum number of rejected steps, pass -1 for unlimited 4613 4614 Notes: 4615 The counter is reset to zero for each step 4616 4617 Options Database Key: 4618 . -ts_max_reject - Maximum number of step rejections before a step fails 4619 4620 Level: intermediate 4621 4622 .keywords: TS, set, maximum, number 4623 4624 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 4625 @*/ 4626 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 4627 { 4628 PetscFunctionBegin; 4629 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4630 ts->max_reject = rejects; 4631 PetscFunctionReturn(0); 4632 } 4633 4634 #undef __FUNCT__ 4635 #define __FUNCT__ "TSSetMaxSNESFailures" 4636 /*@ 4637 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 4638 4639 Not Collective 4640 4641 Input Parameter: 4642 + ts - TS context 4643 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 4644 4645 Notes: 4646 The counter is reset to zero for each successive call to TSSolve(). 4647 4648 Options Database Key: 4649 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 4650 4651 Level: intermediate 4652 4653 .keywords: TS, set, maximum, number 4654 4655 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 4656 @*/ 4657 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 4658 { 4659 PetscFunctionBegin; 4660 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4661 ts->max_snes_failures = fails; 4662 PetscFunctionReturn(0); 4663 } 4664 4665 #undef __FUNCT__ 4666 #define __FUNCT__ "TSSetErrorIfStepFails" 4667 /*@ 4668 TSSetErrorIfStepFails - Error if no step succeeds 4669 4670 Not Collective 4671 4672 Input Parameter: 4673 + ts - TS context 4674 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 4675 4676 Options Database Key: 4677 . -ts_error_if_step_fails - Error if no step succeeds 4678 4679 Level: intermediate 4680 4681 .keywords: TS, set, error 4682 4683 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 4684 @*/ 4685 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 4686 { 4687 PetscFunctionBegin; 4688 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4689 ts->errorifstepfailed = err; 4690 PetscFunctionReturn(0); 4691 } 4692 4693 #undef __FUNCT__ 4694 #define __FUNCT__ "TSMonitorSolutionBinary" 4695 /*@C 4696 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 4697 4698 Collective on TS 4699 4700 Input Parameters: 4701 + ts - the TS context 4702 . step - current time-step 4703 . ptime - current time 4704 . u - current state 4705 - viewer - binary viewer 4706 4707 Level: intermediate 4708 4709 .keywords: TS, vector, monitor, view 4710 4711 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4712 @*/ 4713 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec u,void *viewer) 4714 { 4715 PetscErrorCode ierr; 4716 PetscViewer v = (PetscViewer)viewer; 4717 4718 PetscFunctionBegin; 4719 ierr = VecView(u,v);CHKERRQ(ierr); 4720 PetscFunctionReturn(0); 4721 } 4722 4723 #undef __FUNCT__ 4724 #define __FUNCT__ "TSMonitorSolutionVTK" 4725 /*@C 4726 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 4727 4728 Collective on TS 4729 4730 Input Parameters: 4731 + ts - the TS context 4732 . step - current time-step 4733 . ptime - current time 4734 . u - current state 4735 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4736 4737 Level: intermediate 4738 4739 Notes: 4740 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. 4741 These are named according to the file name template. 4742 4743 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 4744 4745 .keywords: TS, vector, monitor, view 4746 4747 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4748 @*/ 4749 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate) 4750 { 4751 PetscErrorCode ierr; 4752 char filename[PETSC_MAX_PATH_LEN]; 4753 PetscViewer viewer; 4754 4755 PetscFunctionBegin; 4756 ierr = PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);CHKERRQ(ierr); 4757 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 4758 ierr = VecView(u,viewer);CHKERRQ(ierr); 4759 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4760 PetscFunctionReturn(0); 4761 } 4762 4763 #undef __FUNCT__ 4764 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 4765 /*@C 4766 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 4767 4768 Collective on TS 4769 4770 Input Parameters: 4771 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4772 4773 Level: intermediate 4774 4775 Note: 4776 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 4777 4778 .keywords: TS, vector, monitor, view 4779 4780 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 4781 @*/ 4782 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 4783 { 4784 PetscErrorCode ierr; 4785 4786 PetscFunctionBegin; 4787 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 4788 PetscFunctionReturn(0); 4789 } 4790 4791 #undef __FUNCT__ 4792 #define __FUNCT__ "TSGetAdapt" 4793 /*@ 4794 TSGetAdapt - Get the adaptive controller context for the current method 4795 4796 Collective on TS if controller has not been created yet 4797 4798 Input Arguments: 4799 . ts - time stepping context 4800 4801 Output Arguments: 4802 . adapt - adaptive controller 4803 4804 Level: intermediate 4805 4806 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 4807 @*/ 4808 PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt) 4809 { 4810 PetscErrorCode ierr; 4811 4812 PetscFunctionBegin; 4813 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4814 PetscValidPointer(adapt,2); 4815 if (!ts->adapt) { 4816 ierr = TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);CHKERRQ(ierr); 4817 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)ts->adapt);CHKERRQ(ierr); 4818 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 4819 } 4820 *adapt = ts->adapt; 4821 PetscFunctionReturn(0); 4822 } 4823 4824 #undef __FUNCT__ 4825 #define __FUNCT__ "TSSetTolerances" 4826 /*@ 4827 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 4828 4829 Logically Collective 4830 4831 Input Arguments: 4832 + ts - time integration context 4833 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 4834 . vatol - vector of absolute tolerances or NULL, used in preference to atol if present 4835 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 4836 - vrtol - vector of relative tolerances or NULL, used in preference to atol if present 4837 4838 Options Database keys: 4839 + -ts_rtol <rtol> - relative tolerance for local truncation error 4840 - -ts_atol <atol> Absolute tolerance for local truncation error 4841 4842 Level: beginner 4843 4844 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 4845 @*/ 4846 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 4847 { 4848 PetscErrorCode ierr; 4849 4850 PetscFunctionBegin; 4851 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 4852 if (vatol) { 4853 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 4854 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 4855 4856 ts->vatol = vatol; 4857 } 4858 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 4859 if (vrtol) { 4860 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 4861 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 4862 4863 ts->vrtol = vrtol; 4864 } 4865 PetscFunctionReturn(0); 4866 } 4867 4868 #undef __FUNCT__ 4869 #define __FUNCT__ "TSGetTolerances" 4870 /*@ 4871 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 4872 4873 Logically Collective 4874 4875 Input Arguments: 4876 . ts - time integration context 4877 4878 Output Arguments: 4879 + atol - scalar absolute tolerances, NULL to ignore 4880 . vatol - vector of absolute tolerances, NULL to ignore 4881 . rtol - scalar relative tolerances, NULL to ignore 4882 - vrtol - vector of relative tolerances, NULL to ignore 4883 4884 Level: beginner 4885 4886 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 4887 @*/ 4888 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 4889 { 4890 PetscFunctionBegin; 4891 if (atol) *atol = ts->atol; 4892 if (vatol) *vatol = ts->vatol; 4893 if (rtol) *rtol = ts->rtol; 4894 if (vrtol) *vrtol = ts->vrtol; 4895 PetscFunctionReturn(0); 4896 } 4897 4898 #undef __FUNCT__ 4899 #define __FUNCT__ "TSErrorNormWRMS" 4900 /*@ 4901 TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state 4902 4903 Collective on TS 4904 4905 Input Arguments: 4906 + ts - time stepping context 4907 - Y - state vector to be compared to ts->vec_sol 4908 4909 Output Arguments: 4910 . norm - weighted norm, a value of 1.0 is considered small 4911 4912 Level: developer 4913 4914 .seealso: TSSetTolerances() 4915 @*/ 4916 PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm) 4917 { 4918 PetscErrorCode ierr; 4919 PetscInt i,n,N; 4920 const PetscScalar *u,*y; 4921 Vec U; 4922 PetscReal sum,gsum; 4923 4924 PetscFunctionBegin; 4925 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4926 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 4927 PetscValidPointer(norm,3); 4928 U = ts->vec_sol; 4929 PetscCheckSameTypeAndComm(U,1,Y,2); 4930 if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 4931 4932 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 4933 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 4934 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 4935 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 4936 sum = 0.; 4937 if (ts->vatol && ts->vrtol) { 4938 const PetscScalar *atol,*rtol; 4939 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4940 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4941 for (i=0; i<n; i++) { 4942 PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4943 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4944 } 4945 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4946 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4947 } else if (ts->vatol) { /* vector atol, scalar rtol */ 4948 const PetscScalar *atol; 4949 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4950 for (i=0; i<n; i++) { 4951 PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4952 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4953 } 4954 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4955 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 4956 const PetscScalar *rtol; 4957 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4958 for (i=0; i<n; i++) { 4959 PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4960 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4961 } 4962 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4963 } else { /* scalar atol, scalar rtol */ 4964 for (i=0; i<n; i++) { 4965 PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4966 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4967 } 4968 } 4969 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 4970 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 4971 4972 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4973 *norm = PetscSqrtReal(gsum / N); 4974 if (PetscIsInfOrNanReal(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 4975 PetscFunctionReturn(0); 4976 } 4977 4978 #undef __FUNCT__ 4979 #define __FUNCT__ "TSSetCFLTimeLocal" 4980 /*@ 4981 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 4982 4983 Logically Collective on TS 4984 4985 Input Arguments: 4986 + ts - time stepping context 4987 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 4988 4989 Note: 4990 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 4991 4992 Level: intermediate 4993 4994 .seealso: TSGetCFLTime(), TSADAPTCFL 4995 @*/ 4996 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 4997 { 4998 PetscFunctionBegin; 4999 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5000 ts->cfltime_local = cfltime; 5001 ts->cfltime = -1.; 5002 PetscFunctionReturn(0); 5003 } 5004 5005 #undef __FUNCT__ 5006 #define __FUNCT__ "TSGetCFLTime" 5007 /*@ 5008 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 5009 5010 Collective on TS 5011 5012 Input Arguments: 5013 . ts - time stepping context 5014 5015 Output Arguments: 5016 . cfltime - maximum stable time step for forward Euler 5017 5018 Level: advanced 5019 5020 .seealso: TSSetCFLTimeLocal() 5021 @*/ 5022 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 5023 { 5024 PetscErrorCode ierr; 5025 5026 PetscFunctionBegin; 5027 if (ts->cfltime < 0) { 5028 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 5029 } 5030 *cfltime = ts->cfltime; 5031 PetscFunctionReturn(0); 5032 } 5033 5034 #undef __FUNCT__ 5035 #define __FUNCT__ "TSVISetVariableBounds" 5036 /*@ 5037 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 5038 5039 Input Parameters: 5040 . ts - the TS context. 5041 . xl - lower bound. 5042 . xu - upper bound. 5043 5044 Notes: 5045 If this routine is not called then the lower and upper bounds are set to 5046 PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp(). 5047 5048 Level: advanced 5049 5050 @*/ 5051 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 5052 { 5053 PetscErrorCode ierr; 5054 SNES snes; 5055 5056 PetscFunctionBegin; 5057 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 5058 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 5059 PetscFunctionReturn(0); 5060 } 5061 5062 #if defined(PETSC_HAVE_MATLAB_ENGINE) 5063 #include <mex.h> 5064 5065 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 5066 5067 #undef __FUNCT__ 5068 #define __FUNCT__ "TSComputeFunction_Matlab" 5069 /* 5070 TSComputeFunction_Matlab - Calls the function that has been set with 5071 TSSetFunctionMatlab(). 5072 5073 Collective on TS 5074 5075 Input Parameters: 5076 + snes - the TS context 5077 - u - input vector 5078 5079 Output Parameter: 5080 . y - function vector, as set by TSSetFunction() 5081 5082 Notes: 5083 TSComputeFunction() is typically used within nonlinear solvers 5084 implementations, so most users would not generally call this routine 5085 themselves. 5086 5087 Level: developer 5088 5089 .keywords: TS, nonlinear, compute, function 5090 5091 .seealso: TSSetFunction(), TSGetFunction() 5092 */ 5093 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx) 5094 { 5095 PetscErrorCode ierr; 5096 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5097 int nlhs = 1,nrhs = 7; 5098 mxArray *plhs[1],*prhs[7]; 5099 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 5100 5101 PetscFunctionBegin; 5102 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 5103 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 5104 PetscValidHeaderSpecific(udot,VEC_CLASSID,4); 5105 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 5106 PetscCheckSameComm(snes,1,u,3); 5107 PetscCheckSameComm(snes,1,y,5); 5108 5109 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 5110 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5111 ierr = PetscMemcpy(&lxdot,&udot,sizeof(udot));CHKERRQ(ierr); 5112 ierr = PetscMemcpy(&ly,&y,sizeof(u));CHKERRQ(ierr); 5113 5114 prhs[0] = mxCreateDoubleScalar((double)ls); 5115 prhs[1] = mxCreateDoubleScalar(time); 5116 prhs[2] = mxCreateDoubleScalar((double)lx); 5117 prhs[3] = mxCreateDoubleScalar((double)lxdot); 5118 prhs[4] = mxCreateDoubleScalar((double)ly); 5119 prhs[5] = mxCreateString(sctx->funcname); 5120 prhs[6] = sctx->ctx; 5121 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 5122 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5123 mxDestroyArray(prhs[0]); 5124 mxDestroyArray(prhs[1]); 5125 mxDestroyArray(prhs[2]); 5126 mxDestroyArray(prhs[3]); 5127 mxDestroyArray(prhs[4]); 5128 mxDestroyArray(prhs[5]); 5129 mxDestroyArray(plhs[0]); 5130 PetscFunctionReturn(0); 5131 } 5132 5133 5134 #undef __FUNCT__ 5135 #define __FUNCT__ "TSSetFunctionMatlab" 5136 /* 5137 TSSetFunctionMatlab - Sets the function evaluation routine and function 5138 vector for use by the TS routines in solving ODEs 5139 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 5140 5141 Logically Collective on TS 5142 5143 Input Parameters: 5144 + ts - the TS context 5145 - func - function evaluation routine 5146 5147 Calling sequence of func: 5148 $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx); 5149 5150 Level: beginner 5151 5152 .keywords: TS, nonlinear, set, function 5153 5154 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5155 */ 5156 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 5157 { 5158 PetscErrorCode ierr; 5159 TSMatlabContext *sctx; 5160 5161 PetscFunctionBegin; 5162 /* currently sctx is memory bleed */ 5163 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5164 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5165 /* 5166 This should work, but it doesn't 5167 sctx->ctx = ctx; 5168 mexMakeArrayPersistent(sctx->ctx); 5169 */ 5170 sctx->ctx = mxDuplicateArray(ctx); 5171 5172 ierr = TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 5173 PetscFunctionReturn(0); 5174 } 5175 5176 #undef __FUNCT__ 5177 #define __FUNCT__ "TSComputeJacobian_Matlab" 5178 /* 5179 TSComputeJacobian_Matlab - Calls the function that has been set with 5180 TSSetJacobianMatlab(). 5181 5182 Collective on TS 5183 5184 Input Parameters: 5185 + ts - the TS context 5186 . u - input vector 5187 . A, B - the matrices 5188 - ctx - user context 5189 5190 Level: developer 5191 5192 .keywords: TS, nonlinear, compute, function 5193 5194 .seealso: TSSetFunction(), TSGetFunction() 5195 @*/ 5196 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat A,Mat B,void *ctx) 5197 { 5198 PetscErrorCode ierr; 5199 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5200 int nlhs = 2,nrhs = 9; 5201 mxArray *plhs[2],*prhs[9]; 5202 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 5203 5204 PetscFunctionBegin; 5205 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5206 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 5207 5208 /* call Matlab function in ctx with arguments u and y */ 5209 5210 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 5211 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5212 ierr = PetscMemcpy(&lxdot,&udot,sizeof(u));CHKERRQ(ierr); 5213 ierr = PetscMemcpy(&lA,A,sizeof(u));CHKERRQ(ierr); 5214 ierr = PetscMemcpy(&lB,B,sizeof(u));CHKERRQ(ierr); 5215 5216 prhs[0] = mxCreateDoubleScalar((double)ls); 5217 prhs[1] = mxCreateDoubleScalar((double)time); 5218 prhs[2] = mxCreateDoubleScalar((double)lx); 5219 prhs[3] = mxCreateDoubleScalar((double)lxdot); 5220 prhs[4] = mxCreateDoubleScalar((double)shift); 5221 prhs[5] = mxCreateDoubleScalar((double)lA); 5222 prhs[6] = mxCreateDoubleScalar((double)lB); 5223 prhs[7] = mxCreateString(sctx->funcname); 5224 prhs[8] = sctx->ctx; 5225 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 5226 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5227 mxDestroyArray(prhs[0]); 5228 mxDestroyArray(prhs[1]); 5229 mxDestroyArray(prhs[2]); 5230 mxDestroyArray(prhs[3]); 5231 mxDestroyArray(prhs[4]); 5232 mxDestroyArray(prhs[5]); 5233 mxDestroyArray(prhs[6]); 5234 mxDestroyArray(prhs[7]); 5235 mxDestroyArray(plhs[0]); 5236 mxDestroyArray(plhs[1]); 5237 PetscFunctionReturn(0); 5238 } 5239 5240 5241 #undef __FUNCT__ 5242 #define __FUNCT__ "TSSetJacobianMatlab" 5243 /* 5244 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 5245 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 5246 5247 Logically Collective on TS 5248 5249 Input Parameters: 5250 + ts - the TS context 5251 . A,B - Jacobian matrices 5252 . func - function evaluation routine 5253 - ctx - user context 5254 5255 Calling sequence of func: 5256 $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx); 5257 5258 5259 Level: developer 5260 5261 .keywords: TS, nonlinear, set, function 5262 5263 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5264 */ 5265 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 5266 { 5267 PetscErrorCode ierr; 5268 TSMatlabContext *sctx; 5269 5270 PetscFunctionBegin; 5271 /* currently sctx is memory bleed */ 5272 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5273 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5274 /* 5275 This should work, but it doesn't 5276 sctx->ctx = ctx; 5277 mexMakeArrayPersistent(sctx->ctx); 5278 */ 5279 sctx->ctx = mxDuplicateArray(ctx); 5280 5281 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 5282 PetscFunctionReturn(0); 5283 } 5284 5285 #undef __FUNCT__ 5286 #define __FUNCT__ "TSMonitor_Matlab" 5287 /* 5288 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 5289 5290 Collective on TS 5291 5292 .seealso: TSSetFunction(), TSGetFunction() 5293 @*/ 5294 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx) 5295 { 5296 PetscErrorCode ierr; 5297 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 5298 int nlhs = 1,nrhs = 6; 5299 mxArray *plhs[1],*prhs[6]; 5300 long long int lx = 0,ls = 0; 5301 5302 PetscFunctionBegin; 5303 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5304 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 5305 5306 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 5307 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 5308 5309 prhs[0] = mxCreateDoubleScalar((double)ls); 5310 prhs[1] = mxCreateDoubleScalar((double)it); 5311 prhs[2] = mxCreateDoubleScalar((double)time); 5312 prhs[3] = mxCreateDoubleScalar((double)lx); 5313 prhs[4] = mxCreateString(sctx->funcname); 5314 prhs[5] = sctx->ctx; 5315 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 5316 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 5317 mxDestroyArray(prhs[0]); 5318 mxDestroyArray(prhs[1]); 5319 mxDestroyArray(prhs[2]); 5320 mxDestroyArray(prhs[3]); 5321 mxDestroyArray(prhs[4]); 5322 mxDestroyArray(plhs[0]); 5323 PetscFunctionReturn(0); 5324 } 5325 5326 5327 #undef __FUNCT__ 5328 #define __FUNCT__ "TSMonitorSetMatlab" 5329 /* 5330 TSMonitorSetMatlab - Sets the monitor function from Matlab 5331 5332 Level: developer 5333 5334 .keywords: TS, nonlinear, set, function 5335 5336 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 5337 */ 5338 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 5339 { 5340 PetscErrorCode ierr; 5341 TSMatlabContext *sctx; 5342 5343 PetscFunctionBegin; 5344 /* currently sctx is memory bleed */ 5345 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 5346 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 5347 /* 5348 This should work, but it doesn't 5349 sctx->ctx = ctx; 5350 mexMakeArrayPersistent(sctx->ctx); 5351 */ 5352 sctx->ctx = mxDuplicateArray(ctx); 5353 5354 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 5355 PetscFunctionReturn(0); 5356 } 5357 #endif 5358 5359 5360 5361 #undef __FUNCT__ 5362 #define __FUNCT__ "TSMonitorLGSolution" 5363 /*@C 5364 TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector 5365 in a time based line graph 5366 5367 Collective on TS 5368 5369 Input Parameters: 5370 + ts - the TS context 5371 . step - current time-step 5372 . ptime - current time 5373 - lg - a line graph object 5374 5375 Level: intermediate 5376 5377 Notes: each process in a parallel run displays its component solutions in a separate window 5378 5379 .keywords: TS, vector, monitor, view 5380 5381 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 5382 @*/ 5383 PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 5384 { 5385 PetscErrorCode ierr; 5386 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 5387 const PetscScalar *yy; 5388 PetscInt dim; 5389 5390 PetscFunctionBegin; 5391 if (!step) { 5392 PetscDrawAxis axis; 5393 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 5394 ierr = PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");CHKERRQ(ierr); 5395 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 5396 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 5397 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 5398 } 5399 ierr = VecGetArrayRead(u,&yy);CHKERRQ(ierr); 5400 #if defined(PETSC_USE_COMPLEX) 5401 { 5402 PetscReal *yreal; 5403 PetscInt i,n; 5404 ierr = VecGetLocalSize(u,&n);CHKERRQ(ierr); 5405 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 5406 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 5407 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 5408 ierr = PetscFree(yreal);CHKERRQ(ierr); 5409 } 5410 #else 5411 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 5412 #endif 5413 ierr = VecRestoreArrayRead(u,&yy);CHKERRQ(ierr); 5414 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 5415 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 5416 } 5417 PetscFunctionReturn(0); 5418 } 5419 5420 #undef __FUNCT__ 5421 #define __FUNCT__ "TSMonitorLGError" 5422 /*@C 5423 TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector 5424 in a time based line graph 5425 5426 Collective on TS 5427 5428 Input Parameters: 5429 + ts - the TS context 5430 . step - current time-step 5431 . ptime - current time 5432 - lg - a line graph object 5433 5434 Level: intermediate 5435 5436 Notes: 5437 Only for sequential solves. 5438 5439 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 5440 5441 Options Database Keys: 5442 . -ts_monitor_lg_error - create a graphical monitor of error history 5443 5444 .keywords: TS, vector, monitor, view 5445 5446 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 5447 @*/ 5448 PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 5449 { 5450 PetscErrorCode ierr; 5451 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 5452 const PetscScalar *yy; 5453 Vec y; 5454 PetscInt dim; 5455 5456 PetscFunctionBegin; 5457 if (!step) { 5458 PetscDrawAxis axis; 5459 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 5460 ierr = PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");CHKERRQ(ierr); 5461 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 5462 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 5463 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 5464 } 5465 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 5466 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 5467 ierr = VecAXPY(y,-1.0,u);CHKERRQ(ierr); 5468 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 5469 #if defined(PETSC_USE_COMPLEX) 5470 { 5471 PetscReal *yreal; 5472 PetscInt i,n; 5473 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 5474 ierr = PetscMalloc1(n,&yreal);CHKERRQ(ierr); 5475 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 5476 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 5477 ierr = PetscFree(yreal);CHKERRQ(ierr); 5478 } 5479 #else 5480 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 5481 #endif 5482 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 5483 ierr = VecDestroy(&y);CHKERRQ(ierr); 5484 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 5485 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 5486 } 5487 PetscFunctionReturn(0); 5488 } 5489 5490 #undef __FUNCT__ 5491 #define __FUNCT__ "TSMonitorLGSNESIterations" 5492 PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 5493 { 5494 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 5495 PetscReal x = ptime,y; 5496 PetscErrorCode ierr; 5497 PetscInt its; 5498 5499 PetscFunctionBegin; 5500 if (!n) { 5501 PetscDrawAxis axis; 5502 5503 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 5504 ierr = PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");CHKERRQ(ierr); 5505 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 5506 5507 ctx->snes_its = 0; 5508 } 5509 ierr = TSGetSNESIterations(ts,&its);CHKERRQ(ierr); 5510 y = its - ctx->snes_its; 5511 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 5512 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 5513 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 5514 } 5515 ctx->snes_its = its; 5516 PetscFunctionReturn(0); 5517 } 5518 5519 #undef __FUNCT__ 5520 #define __FUNCT__ "TSMonitorLGKSPIterations" 5521 PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 5522 { 5523 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 5524 PetscReal x = ptime,y; 5525 PetscErrorCode ierr; 5526 PetscInt its; 5527 5528 PetscFunctionBegin; 5529 if (!n) { 5530 PetscDrawAxis axis; 5531 5532 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 5533 ierr = PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");CHKERRQ(ierr); 5534 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 5535 5536 ctx->ksp_its = 0; 5537 } 5538 ierr = TSGetKSPIterations(ts,&its);CHKERRQ(ierr); 5539 y = its - ctx->ksp_its; 5540 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 5541 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 5542 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 5543 } 5544 ctx->ksp_its = its; 5545 PetscFunctionReturn(0); 5546 } 5547 5548 #undef __FUNCT__ 5549 #define __FUNCT__ "TSComputeLinearStability" 5550 /*@ 5551 TSComputeLinearStability - computes the linear stability function at a point 5552 5553 Collective on TS and Vec 5554 5555 Input Parameters: 5556 + ts - the TS context 5557 - xr,xi - real and imaginary part of input arguments 5558 5559 Output Parameters: 5560 . yr,yi - real and imaginary part of function value 5561 5562 Level: developer 5563 5564 .keywords: TS, compute 5565 5566 .seealso: TSSetRHSFunction(), TSComputeIFunction() 5567 @*/ 5568 PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 5569 { 5570 PetscErrorCode ierr; 5571 5572 PetscFunctionBegin; 5573 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 5574 if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method"); 5575 ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr); 5576 PetscFunctionReturn(0); 5577 } 5578 5579 #undef __FUNCT__ 5580 #define __FUNCT__ "TSRollBack" 5581 /*@ 5582 TSRollBack - Rolls back one time step 5583 5584 Collective on TS 5585 5586 Input Parameter: 5587 . ts - the TS context obtained from TSCreate() 5588 5589 Level: advanced 5590 5591 .keywords: TS, timestep, rollback 5592 5593 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate() 5594 @*/ 5595 PetscErrorCode TSRollBack(TS ts) 5596 { 5597 PetscErrorCode ierr; 5598 5599 PetscFunctionBegin; 5600 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 5601 5602 if (!ts->ops->rollback) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRollBack not implemented for type '%s'",((PetscObject)ts)->type_name); 5603 ierr = (*ts->ops->rollback)(ts);CHKERRQ(ierr); 5604 ts->time_step = ts->ptime - ts->ptime_prev; 5605 ts->ptime = ts->ptime_prev; 5606 PetscFunctionReturn(0); 5607 } 5608 5609 #undef __FUNCT__ 5610 #define __FUNCT__ "TSGetStages" 5611 /*@ 5612 TSGetStages - Get the number of stages and stage values 5613 5614 Input Parameter: 5615 . ts - the TS context obtained from TSCreate() 5616 5617 Level: advanced 5618 5619 .keywords: TS, getstages 5620 5621 .seealso: TSCreate() 5622 @*/ 5623 PetscErrorCode TSGetStages(TS ts,PetscInt *ns, Vec **Y) 5624 { 5625 PetscErrorCode ierr; 5626 5627 PetscFunctionBegin; 5628 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 5629 PetscValidPointer(ns,2); 5630 5631 if (!ts->ops->getstages) *ns=0; 5632 else { 5633 ierr = (*ts->ops->getstages)(ts,ns,Y);CHKERRQ(ierr); 5634 } 5635 PetscFunctionReturn(0); 5636 } 5637 5638