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