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