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