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 }; 2859 2860 #undef __FUNCT__ 2861 #define __FUNCT__ "TSMonitorDrawSolution" 2862 /*@C 2863 TSMonitorDrawSolution - Monitors progress of the TS solvers by calling 2864 VecView() for the solution at each timestep 2865 2866 Collective on TS 2867 2868 Input Parameters: 2869 + ts - the TS context 2870 . step - current time-step 2871 . ptime - current time 2872 - dummy - either a viewer or PETSC_NULL 2873 2874 Options Database: 2875 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 2876 2877 Notes: the initial solution and current solution are not displayed with a common axis scaling so generally the option -ts_monitor_draw_solution_initial 2878 will look bad 2879 2880 Level: intermediate 2881 2882 .keywords: TS, vector, monitor, view 2883 2884 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 2885 @*/ 2886 PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 2887 { 2888 PetscErrorCode ierr; 2889 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 2890 2891 PetscFunctionBegin; 2892 if (!step && ictx->showinitial) { 2893 if (!ictx->initialsolution) { 2894 ierr = VecDuplicate(u,&ictx->initialsolution);CHKERRQ(ierr); 2895 } 2896 ierr = VecCopy(u,ictx->initialsolution);CHKERRQ(ierr); 2897 } 2898 if (!(((ictx->howoften > 0) && (!(step % ictx->howoften)) && (step > -1)) || ((ictx->howoften == -1) && (step == -1)))) PetscFunctionReturn(0); 2899 2900 if (ictx->showinitial) { 2901 PetscReal pause; 2902 ierr = PetscViewerDrawGetPause(ictx->viewer,&pause);CHKERRQ(ierr); 2903 ierr = PetscViewerDrawSetPause(ictx->viewer,0.0);CHKERRQ(ierr); 2904 ierr = VecView(ictx->initialsolution,ictx->viewer);CHKERRQ(ierr); 2905 ierr = PetscViewerDrawSetPause(ictx->viewer,pause);CHKERRQ(ierr); 2906 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);CHKERRQ(ierr); 2907 } 2908 ierr = VecView(u,ictx->viewer);CHKERRQ(ierr); 2909 if (ictx->showinitial) { 2910 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);CHKERRQ(ierr); 2911 } 2912 PetscFunctionReturn(0); 2913 } 2914 2915 2916 #undef __FUNCT__ 2917 #define __FUNCT__ "TSMonitorDrawCtxDestroy" 2918 /*@C 2919 TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution() 2920 2921 Collective on TS 2922 2923 Input Parameters: 2924 . ctx - the monitor context 2925 2926 Level: intermediate 2927 2928 .keywords: TS, vector, monitor, view 2929 2930 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError() 2931 @*/ 2932 PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx) 2933 { 2934 PetscErrorCode ierr; 2935 2936 PetscFunctionBegin; 2937 ierr = PetscViewerDestroy(&(*ictx)->viewer);CHKERRQ(ierr); 2938 ierr = VecDestroy(&(*ictx)->initialsolution);CHKERRQ(ierr); 2939 ierr = PetscFree(*ictx);CHKERRQ(ierr); 2940 PetscFunctionReturn(0); 2941 } 2942 2943 #undef __FUNCT__ 2944 #define __FUNCT__ "TSMonitorDrawCtxCreate" 2945 /*@C 2946 TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx 2947 2948 Collective on TS 2949 2950 Input Parameter: 2951 . ts - time-step context 2952 2953 Output Patameter: 2954 . ctx - the monitor context 2955 2956 Options Database: 2957 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 2958 2959 Level: intermediate 2960 2961 .keywords: TS, vector, monitor, view 2962 2963 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx() 2964 @*/ 2965 PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx) 2966 { 2967 PetscErrorCode ierr; 2968 2969 PetscFunctionBegin; 2970 ierr = PetscNew(struct _n_TSMonitorDrawCtx,ctx);CHKERRQ(ierr); 2971 ierr = PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);CHKERRQ(ierr); 2972 2973 (*ctx)->showinitial = PETSC_FALSE; 2974 (*ctx)->howoften = howoften; 2975 2976 ierr = PetscOptionsGetBool(PETSC_NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,PETSC_NULL);CHKERRQ(ierr); 2977 PetscFunctionReturn(0); 2978 } 2979 2980 #undef __FUNCT__ 2981 #define __FUNCT__ "TSMonitorDrawError" 2982 /*@C 2983 TSMonitorDrawError - Monitors progress of the TS solvers by calling 2984 VecView() for the error at each timestep 2985 2986 Collective on TS 2987 2988 Input Parameters: 2989 + ts - the TS context 2990 . step - current time-step 2991 . ptime - current time 2992 - dummy - either a viewer or PETSC_NULL 2993 2994 Level: intermediate 2995 2996 .keywords: TS, vector, monitor, view 2997 2998 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 2999 @*/ 3000 PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3001 { 3002 PetscErrorCode ierr; 3003 TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy; 3004 PetscViewer viewer = ctx->viewer; 3005 Vec work; 3006 3007 PetscFunctionBegin; 3008 if (!(((ctx->howoften > 0) && (!(step % ctx->howoften)) && (step > -1)) || ((ctx->howoften == -1) && (step == -1)))) PetscFunctionReturn(0); 3009 ierr = VecDuplicate(u,&work);CHKERRQ(ierr); 3010 ierr = TSComputeSolutionFunction(ts,ptime,work);CHKERRQ(ierr); 3011 ierr = VecAXPY(work,-1.0,u);CHKERRQ(ierr); 3012 ierr = VecView(work,viewer);CHKERRQ(ierr); 3013 ierr = VecDestroy(&work);CHKERRQ(ierr); 3014 PetscFunctionReturn(0); 3015 } 3016 3017 #include <petsc-private/dmimpl.h> 3018 #undef __FUNCT__ 3019 #define __FUNCT__ "TSSetDM" 3020 /*@ 3021 TSSetDM - Sets the DM that may be used by some preconditioners 3022 3023 Logically Collective on TS and DM 3024 3025 Input Parameters: 3026 + ts - the preconditioner context 3027 - dm - the dm 3028 3029 Level: intermediate 3030 3031 3032 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 3033 @*/ 3034 PetscErrorCode TSSetDM(TS ts,DM dm) 3035 { 3036 PetscErrorCode ierr; 3037 SNES snes; 3038 DMTS tsdm; 3039 3040 PetscFunctionBegin; 3041 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3042 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 3043 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 3044 if (ts->dm->dmts && !dm->dmts) { 3045 ierr = DMCopyDMTS(ts->dm,dm);CHKERRQ(ierr); 3046 ierr = DMGetDMTS(ts->dm,&tsdm);CHKERRQ(ierr); 3047 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 3048 tsdm->originaldm = dm; 3049 } 3050 } 3051 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 3052 } 3053 ts->dm = dm; 3054 3055 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3056 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 3057 PetscFunctionReturn(0); 3058 } 3059 3060 #undef __FUNCT__ 3061 #define __FUNCT__ "TSGetDM" 3062 /*@ 3063 TSGetDM - Gets the DM that may be used by some preconditioners 3064 3065 Not Collective 3066 3067 Input Parameter: 3068 . ts - the preconditioner context 3069 3070 Output Parameter: 3071 . dm - the dm 3072 3073 Level: intermediate 3074 3075 3076 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 3077 @*/ 3078 PetscErrorCode TSGetDM(TS ts,DM *dm) 3079 { 3080 PetscErrorCode ierr; 3081 3082 PetscFunctionBegin; 3083 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3084 if (!ts->dm) { 3085 ierr = DMShellCreate(((PetscObject)ts)->comm,&ts->dm);CHKERRQ(ierr); 3086 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 3087 } 3088 *dm = ts->dm; 3089 PetscFunctionReturn(0); 3090 } 3091 3092 #undef __FUNCT__ 3093 #define __FUNCT__ "SNESTSFormFunction" 3094 /*@ 3095 SNESTSFormFunction - Function to evaluate nonlinear residual 3096 3097 Logically Collective on SNES 3098 3099 Input Parameter: 3100 + snes - nonlinear solver 3101 . U - the current state at which to evaluate the residual 3102 - ctx - user context, must be a TS 3103 3104 Output Parameter: 3105 . F - the nonlinear residual 3106 3107 Notes: 3108 This function is not normally called by users and is automatically registered with the SNES used by TS. 3109 It is most frequently passed to MatFDColoringSetFunction(). 3110 3111 Level: advanced 3112 3113 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 3114 @*/ 3115 PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx) 3116 { 3117 TS ts = (TS)ctx; 3118 PetscErrorCode ierr; 3119 3120 PetscFunctionBegin; 3121 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3122 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3123 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 3124 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 3125 ierr = (ts->ops->snesfunction)(snes,U,F,ts);CHKERRQ(ierr); 3126 PetscFunctionReturn(0); 3127 } 3128 3129 #undef __FUNCT__ 3130 #define __FUNCT__ "SNESTSFormJacobian" 3131 /*@ 3132 SNESTSFormJacobian - Function to evaluate the Jacobian 3133 3134 Collective on SNES 3135 3136 Input Parameter: 3137 + snes - nonlinear solver 3138 . U - the current state at which to evaluate the residual 3139 - ctx - user context, must be a TS 3140 3141 Output Parameter: 3142 + A - the Jacobian 3143 . B - the preconditioning matrix (may be the same as A) 3144 - flag - indicates any structure change in the matrix 3145 3146 Notes: 3147 This function is not normally called by users and is automatically registered with the SNES used by TS. 3148 3149 Level: developer 3150 3151 .seealso: SNESSetJacobian() 3152 @*/ 3153 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat *A,Mat *B,MatStructure *flag,void *ctx) 3154 { 3155 TS ts = (TS)ctx; 3156 PetscErrorCode ierr; 3157 3158 PetscFunctionBegin; 3159 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3160 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3161 PetscValidPointer(A,3); 3162 PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 3163 PetscValidPointer(B,4); 3164 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); 3165 PetscValidPointer(flag,5); 3166 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 3167 ierr = (ts->ops->snesjacobian)(snes,U,A,B,flag,ts);CHKERRQ(ierr); 3168 PetscFunctionReturn(0); 3169 } 3170 3171 #undef __FUNCT__ 3172 #define __FUNCT__ "TSComputeRHSFunctionLinear" 3173 /*@C 3174 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 3175 3176 Collective on TS 3177 3178 Input Arguments: 3179 + ts - time stepping context 3180 . t - time at which to evaluate 3181 . U - state at which to evaluate 3182 - ctx - context 3183 3184 Output Arguments: 3185 . F - right hand side 3186 3187 Level: intermediate 3188 3189 Notes: 3190 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 3191 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 3192 3193 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 3194 @*/ 3195 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 3196 { 3197 PetscErrorCode ierr; 3198 Mat Arhs,Brhs; 3199 MatStructure flg2; 3200 3201 PetscFunctionBegin; 3202 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 3203 ierr = TSComputeRHSJacobian(ts,t,U,&Arhs,&Brhs,&flg2);CHKERRQ(ierr); 3204 ierr = MatMult(Arhs,U,F);CHKERRQ(ierr); 3205 PetscFunctionReturn(0); 3206 } 3207 3208 #undef __FUNCT__ 3209 #define __FUNCT__ "TSComputeRHSJacobianConstant" 3210 /*@C 3211 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 3212 3213 Collective on TS 3214 3215 Input Arguments: 3216 + ts - time stepping context 3217 . t - time at which to evaluate 3218 . U - state at which to evaluate 3219 - ctx - context 3220 3221 Output Arguments: 3222 + A - pointer to operator 3223 . B - pointer to preconditioning matrix 3224 - flg - matrix structure flag 3225 3226 Level: intermediate 3227 3228 Notes: 3229 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 3230 3231 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 3232 @*/ 3233 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat *A,Mat *B,MatStructure *flg,void *ctx) 3234 { 3235 PetscFunctionBegin; 3236 *flg = SAME_PRECONDITIONER; 3237 PetscFunctionReturn(0); 3238 } 3239 3240 #undef __FUNCT__ 3241 #define __FUNCT__ "TSComputeIFunctionLinear" 3242 /*@C 3243 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 3244 3245 Collective on TS 3246 3247 Input Arguments: 3248 + ts - time stepping context 3249 . t - time at which to evaluate 3250 . U - state at which to evaluate 3251 . Udot - time derivative of state vector 3252 - ctx - context 3253 3254 Output Arguments: 3255 . F - left hand side 3256 3257 Level: intermediate 3258 3259 Notes: 3260 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 3261 user is required to write their own TSComputeIFunction. 3262 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 3263 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 3264 3265 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 3266 @*/ 3267 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 3268 { 3269 PetscErrorCode ierr; 3270 Mat A,B; 3271 MatStructure flg2; 3272 3273 PetscFunctionBegin; 3274 ierr = TSGetIJacobian(ts,&A,&B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 3275 ierr = TSComputeIJacobian(ts,t,U,Udot,1.0,&A,&B,&flg2,PETSC_TRUE);CHKERRQ(ierr); 3276 ierr = MatMult(A,Udot,F);CHKERRQ(ierr); 3277 PetscFunctionReturn(0); 3278 } 3279 3280 #undef __FUNCT__ 3281 #define __FUNCT__ "TSComputeIJacobianConstant" 3282 /*@C 3283 TSComputeIJacobianConstant - Reuses a Jacobian that is time-independent. 3284 3285 Collective on TS 3286 3287 Input Arguments: 3288 + ts - time stepping context 3289 . t - time at which to evaluate 3290 . U - state at which to evaluate 3291 . Udot - time derivative of state vector 3292 . shift - shift to apply 3293 - ctx - context 3294 3295 Output Arguments: 3296 + A - pointer to operator 3297 . B - pointer to preconditioning matrix 3298 - flg - matrix structure flag 3299 3300 Level: intermediate 3301 3302 Notes: 3303 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 3304 3305 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 3306 @*/ 3307 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,void *ctx) 3308 { 3309 PetscFunctionBegin; 3310 *flg = SAME_PRECONDITIONER; 3311 PetscFunctionReturn(0); 3312 } 3313 #undef __FUNCT__ 3314 #define __FUNCT__ "TSGetEquationType" 3315 /*@ 3316 TSGetEquationType - Gets the type of the equation that TS is solving. 3317 3318 Not Collective 3319 3320 Input Parameter: 3321 . ts - the TS context 3322 3323 Output Parameter: 3324 . equation_type - see TSEquatioType 3325 3326 Level: beginner 3327 3328 .keywords: TS, equation type 3329 3330 .seealso: TSSetEquationType(), TSEquationType 3331 @*/ 3332 PetscErrorCode TSGetEquationType(TS ts,TSEquationType *equation_type) 3333 { 3334 PetscFunctionBegin; 3335 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3336 PetscValidPointer(equation_type,2); 3337 *equation_type = ts->equation_type; 3338 PetscFunctionReturn(0); 3339 } 3340 3341 #undef __FUNCT__ 3342 #define __FUNCT__ "TSSetEquationType" 3343 /*@ 3344 TSSetEquationType - Sets the type of the equation that TS is solving. 3345 3346 Not Collective 3347 3348 Input Parameter: 3349 + ts - the TS context 3350 . equation_type - see TSEquatioType 3351 3352 Level: advanced 3353 3354 .keywords: TS, equation type 3355 3356 .seealso: TSGetEquationType(), TSEquationType 3357 @*/ 3358 PetscErrorCode TSSetEquationType(TS ts,TSEquationType equation_type) 3359 { 3360 PetscFunctionBegin; 3361 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3362 ts->equation_type = equation_type; 3363 PetscFunctionReturn(0); 3364 } 3365 3366 #undef __FUNCT__ 3367 #define __FUNCT__ "TSGetConvergedReason" 3368 /*@ 3369 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 3370 3371 Not Collective 3372 3373 Input Parameter: 3374 . ts - the TS context 3375 3376 Output Parameter: 3377 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3378 manual pages for the individual convergence tests for complete lists 3379 3380 Level: beginner 3381 3382 Notes: 3383 Can only be called after the call to TSSolve() is complete. 3384 3385 .keywords: TS, nonlinear, set, convergence, test 3386 3387 .seealso: TSSetConvergenceTest(), TSConvergedReason 3388 @*/ 3389 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 3390 { 3391 PetscFunctionBegin; 3392 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3393 PetscValidPointer(reason,2); 3394 *reason = ts->reason; 3395 PetscFunctionReturn(0); 3396 } 3397 3398 #undef __FUNCT__ 3399 #define __FUNCT__ "TSSetConvergedReason" 3400 /*@ 3401 TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve. 3402 3403 Not Collective 3404 3405 Input Parameter: 3406 + ts - the TS context 3407 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3408 manual pages for the individual convergence tests for complete lists 3409 3410 Level: advanced 3411 3412 Notes: 3413 Can only be called during TSSolve() is active. 3414 3415 .keywords: TS, nonlinear, set, convergence, test 3416 3417 .seealso: TSConvergedReason 3418 @*/ 3419 PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason) 3420 { 3421 PetscFunctionBegin; 3422 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3423 ts->reason = reason; 3424 PetscFunctionReturn(0); 3425 } 3426 3427 #undef __FUNCT__ 3428 #define __FUNCT__ "TSGetSolveTime" 3429 /*@ 3430 TSGetSolveTime - Gets the time after a call to TSSolve() 3431 3432 Not Collective 3433 3434 Input Parameter: 3435 . ts - the TS context 3436 3437 Output Parameter: 3438 . ftime - the final time. This time should correspond to the final time set with TSSetDuration() 3439 3440 Level: beginner 3441 3442 Notes: 3443 Can only be called after the call to TSSolve() is complete. 3444 3445 .keywords: TS, nonlinear, set, convergence, test 3446 3447 .seealso: TSSetConvergenceTest(), TSConvergedReason 3448 @*/ 3449 PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime) 3450 { 3451 PetscFunctionBegin; 3452 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3453 PetscValidPointer(ftime,2); 3454 *ftime = ts->solvetime; 3455 PetscFunctionReturn(0); 3456 } 3457 3458 #undef __FUNCT__ 3459 #define __FUNCT__ "TSGetSNESIterations" 3460 /*@ 3461 TSGetSNESIterations - Gets the total number of nonlinear iterations 3462 used by the time integrator. 3463 3464 Not Collective 3465 3466 Input Parameter: 3467 . ts - TS context 3468 3469 Output Parameter: 3470 . nits - number of nonlinear iterations 3471 3472 Notes: 3473 This counter is reset to zero for each successive call to TSSolve(). 3474 3475 Level: intermediate 3476 3477 .keywords: TS, get, number, nonlinear, iterations 3478 3479 .seealso: TSGetKSPIterations() 3480 @*/ 3481 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 3482 { 3483 PetscFunctionBegin; 3484 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3485 PetscValidIntPointer(nits,2); 3486 *nits = ts->snes_its; 3487 PetscFunctionReturn(0); 3488 } 3489 3490 #undef __FUNCT__ 3491 #define __FUNCT__ "TSGetKSPIterations" 3492 /*@ 3493 TSGetKSPIterations - Gets the total number of linear iterations 3494 used by the time integrator. 3495 3496 Not Collective 3497 3498 Input Parameter: 3499 . ts - TS context 3500 3501 Output Parameter: 3502 . lits - number of linear iterations 3503 3504 Notes: 3505 This counter is reset to zero for each successive call to TSSolve(). 3506 3507 Level: intermediate 3508 3509 .keywords: TS, get, number, linear, iterations 3510 3511 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 3512 @*/ 3513 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 3514 { 3515 PetscFunctionBegin; 3516 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3517 PetscValidIntPointer(lits,2); 3518 *lits = ts->ksp_its; 3519 PetscFunctionReturn(0); 3520 } 3521 3522 #undef __FUNCT__ 3523 #define __FUNCT__ "TSGetStepRejections" 3524 /*@ 3525 TSGetStepRejections - Gets the total number of rejected steps. 3526 3527 Not Collective 3528 3529 Input Parameter: 3530 . ts - TS context 3531 3532 Output Parameter: 3533 . rejects - number of steps rejected 3534 3535 Notes: 3536 This counter is reset to zero for each successive call to TSSolve(). 3537 3538 Level: intermediate 3539 3540 .keywords: TS, get, number 3541 3542 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 3543 @*/ 3544 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 3545 { 3546 PetscFunctionBegin; 3547 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3548 PetscValidIntPointer(rejects,2); 3549 *rejects = ts->reject; 3550 PetscFunctionReturn(0); 3551 } 3552 3553 #undef __FUNCT__ 3554 #define __FUNCT__ "TSGetSNESFailures" 3555 /*@ 3556 TSGetSNESFailures - Gets the total number of failed SNES solves 3557 3558 Not Collective 3559 3560 Input Parameter: 3561 . ts - TS context 3562 3563 Output Parameter: 3564 . fails - number of failed nonlinear solves 3565 3566 Notes: 3567 This counter is reset to zero for each successive call to TSSolve(). 3568 3569 Level: intermediate 3570 3571 .keywords: TS, get, number 3572 3573 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 3574 @*/ 3575 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 3576 { 3577 PetscFunctionBegin; 3578 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3579 PetscValidIntPointer(fails,2); 3580 *fails = ts->num_snes_failures; 3581 PetscFunctionReturn(0); 3582 } 3583 3584 #undef __FUNCT__ 3585 #define __FUNCT__ "TSSetMaxStepRejections" 3586 /*@ 3587 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 3588 3589 Not Collective 3590 3591 Input Parameter: 3592 + ts - TS context 3593 - rejects - maximum number of rejected steps, pass -1 for unlimited 3594 3595 Notes: 3596 The counter is reset to zero for each step 3597 3598 Options Database Key: 3599 . -ts_max_reject - Maximum number of step rejections before a step fails 3600 3601 Level: intermediate 3602 3603 .keywords: TS, set, maximum, number 3604 3605 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 3606 @*/ 3607 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 3608 { 3609 PetscFunctionBegin; 3610 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3611 ts->max_reject = rejects; 3612 PetscFunctionReturn(0); 3613 } 3614 3615 #undef __FUNCT__ 3616 #define __FUNCT__ "TSSetMaxSNESFailures" 3617 /*@ 3618 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 3619 3620 Not Collective 3621 3622 Input Parameter: 3623 + ts - TS context 3624 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 3625 3626 Notes: 3627 The counter is reset to zero for each successive call to TSSolve(). 3628 3629 Options Database Key: 3630 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 3631 3632 Level: intermediate 3633 3634 .keywords: TS, set, maximum, number 3635 3636 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 3637 @*/ 3638 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 3639 { 3640 PetscFunctionBegin; 3641 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3642 ts->max_snes_failures = fails; 3643 PetscFunctionReturn(0); 3644 } 3645 3646 #undef __FUNCT__ 3647 #define __FUNCT__ "TSSetErrorIfStepFails()" 3648 /*@ 3649 TSSetErrorIfStepFails - Error if no step succeeds 3650 3651 Not Collective 3652 3653 Input Parameter: 3654 + ts - TS context 3655 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 3656 3657 Options Database Key: 3658 . -ts_error_if_step_fails - Error if no step succeeds 3659 3660 Level: intermediate 3661 3662 .keywords: TS, set, error 3663 3664 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 3665 @*/ 3666 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 3667 { 3668 PetscFunctionBegin; 3669 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3670 ts->errorifstepfailed = err; 3671 PetscFunctionReturn(0); 3672 } 3673 3674 #undef __FUNCT__ 3675 #define __FUNCT__ "TSMonitorSolutionBinary" 3676 /*@C 3677 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 3678 3679 Collective on TS 3680 3681 Input Parameters: 3682 + ts - the TS context 3683 . step - current time-step 3684 . ptime - current time 3685 . u - current state 3686 - viewer - binary viewer 3687 3688 Level: intermediate 3689 3690 .keywords: TS, vector, monitor, view 3691 3692 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3693 @*/ 3694 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec u,void *viewer) 3695 { 3696 PetscErrorCode ierr; 3697 PetscViewer v = (PetscViewer)viewer; 3698 3699 PetscFunctionBegin; 3700 ierr = VecView(u,v);CHKERRQ(ierr); 3701 PetscFunctionReturn(0); 3702 } 3703 3704 #undef __FUNCT__ 3705 #define __FUNCT__ "TSMonitorSolutionVTK" 3706 /*@C 3707 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 3708 3709 Collective on TS 3710 3711 Input Parameters: 3712 + ts - the TS context 3713 . step - current time-step 3714 . ptime - current time 3715 . u - current state 3716 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 3717 3718 Level: intermediate 3719 3720 Notes: 3721 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. 3722 These are named according to the file name template. 3723 3724 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 3725 3726 .keywords: TS, vector, monitor, view 3727 3728 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3729 @*/ 3730 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate) 3731 { 3732 PetscErrorCode ierr; 3733 char filename[PETSC_MAX_PATH_LEN]; 3734 PetscViewer viewer; 3735 3736 PetscFunctionBegin; 3737 ierr = PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);CHKERRQ(ierr); 3738 ierr = PetscViewerVTKOpen(((PetscObject)ts)->comm,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 3739 ierr = VecView(u,viewer);CHKERRQ(ierr); 3740 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3741 PetscFunctionReturn(0); 3742 } 3743 3744 #undef __FUNCT__ 3745 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 3746 /*@C 3747 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 3748 3749 Collective on TS 3750 3751 Input Parameters: 3752 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 3753 3754 Level: intermediate 3755 3756 Note: 3757 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 3758 3759 .keywords: TS, vector, monitor, view 3760 3761 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 3762 @*/ 3763 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 3764 { 3765 PetscErrorCode ierr; 3766 3767 PetscFunctionBegin; 3768 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 3769 PetscFunctionReturn(0); 3770 } 3771 3772 #undef __FUNCT__ 3773 #define __FUNCT__ "TSGetTSAdapt" 3774 /*@ 3775 TSGetTSAdapt - Get the adaptive controller context for the current method 3776 3777 Collective on TS if controller has not been created yet 3778 3779 Input Arguments: 3780 . ts - time stepping context 3781 3782 Output Arguments: 3783 . adapt - adaptive controller 3784 3785 Level: intermediate 3786 3787 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 3788 @*/ 3789 PetscErrorCode TSGetTSAdapt(TS ts,TSAdapt *adapt) 3790 { 3791 PetscErrorCode ierr; 3792 3793 PetscFunctionBegin; 3794 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3795 PetscValidPointer(adapt,2); 3796 if (!ts->adapt) { 3797 ierr = TSAdaptCreate(((PetscObject)ts)->comm,&ts->adapt);CHKERRQ(ierr); 3798 ierr = PetscLogObjectParent(ts,ts->adapt);CHKERRQ(ierr); 3799 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 3800 } 3801 *adapt = ts->adapt; 3802 PetscFunctionReturn(0); 3803 } 3804 3805 #undef __FUNCT__ 3806 #define __FUNCT__ "TSSetTolerances" 3807 /*@ 3808 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 3809 3810 Logically Collective 3811 3812 Input Arguments: 3813 + ts - time integration context 3814 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 3815 . vatol - vector of absolute tolerances or PETSC_NULL, used in preference to atol if present 3816 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 3817 - vrtol - vector of relative tolerances or PETSC_NULL, used in preference to atol if present 3818 3819 Level: beginner 3820 3821 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 3822 @*/ 3823 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 3824 { 3825 PetscErrorCode ierr; 3826 3827 PetscFunctionBegin; 3828 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 3829 if (vatol) { 3830 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 3831 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 3832 3833 ts->vatol = vatol; 3834 } 3835 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 3836 if (vrtol) { 3837 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 3838 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 3839 3840 ts->vrtol = vrtol; 3841 } 3842 PetscFunctionReturn(0); 3843 } 3844 3845 #undef __FUNCT__ 3846 #define __FUNCT__ "TSGetTolerances" 3847 /*@ 3848 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 3849 3850 Logically Collective 3851 3852 Input Arguments: 3853 . ts - time integration context 3854 3855 Output Arguments: 3856 + atol - scalar absolute tolerances, PETSC_NULL to ignore 3857 . vatol - vector of absolute tolerances, PETSC_NULL to ignore 3858 . rtol - scalar relative tolerances, PETSC_NULL to ignore 3859 - vrtol - vector of relative tolerances, PETSC_NULL to ignore 3860 3861 Level: beginner 3862 3863 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 3864 @*/ 3865 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 3866 { 3867 PetscFunctionBegin; 3868 if (atol) *atol = ts->atol; 3869 if (vatol) *vatol = ts->vatol; 3870 if (rtol) *rtol = ts->rtol; 3871 if (vrtol) *vrtol = ts->vrtol; 3872 PetscFunctionReturn(0); 3873 } 3874 3875 #undef __FUNCT__ 3876 #define __FUNCT__ "TSErrorNormWRMS" 3877 /*@ 3878 TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state 3879 3880 Collective on TS 3881 3882 Input Arguments: 3883 + ts - time stepping context 3884 - Y - state vector to be compared to ts->vec_sol 3885 3886 Output Arguments: 3887 . norm - weighted norm, a value of 1.0 is considered small 3888 3889 Level: developer 3890 3891 .seealso: TSSetTolerances() 3892 @*/ 3893 PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm) 3894 { 3895 PetscErrorCode ierr; 3896 PetscInt i,n,N; 3897 const PetscScalar *u,*y; 3898 Vec U; 3899 PetscReal sum,gsum; 3900 3901 PetscFunctionBegin; 3902 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3903 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 3904 PetscValidPointer(norm,3); 3905 U = ts->vec_sol; 3906 PetscCheckSameTypeAndComm(U,1,Y,2); 3907 if (U == Y) SETERRQ(((PetscObject)U)->comm,PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 3908 3909 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 3910 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 3911 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 3912 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 3913 sum = 0.; 3914 if (ts->vatol && ts->vrtol) { 3915 const PetscScalar *atol,*rtol; 3916 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3917 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3918 for (i=0; i<n; i++) { 3919 PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3920 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3921 } 3922 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3923 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3924 } else if (ts->vatol) { /* vector atol, scalar rtol */ 3925 const PetscScalar *atol; 3926 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3927 for (i=0; i<n; i++) { 3928 PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3929 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3930 } 3931 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 3932 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 3933 const PetscScalar *rtol; 3934 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3935 for (i=0; i<n; i++) { 3936 PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3937 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3938 } 3939 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 3940 } else { /* scalar atol, scalar rtol */ 3941 for (i=0; i<n; i++) { 3942 PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 3943 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 3944 } 3945 } 3946 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 3947 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 3948 3949 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,((PetscObject)ts)->comm);CHKERRQ(ierr); 3950 *norm = PetscSqrtReal(gsum / N); 3951 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 3952 PetscFunctionReturn(0); 3953 } 3954 3955 #undef __FUNCT__ 3956 #define __FUNCT__ "TSSetCFLTimeLocal" 3957 /*@ 3958 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 3959 3960 Logically Collective on TS 3961 3962 Input Arguments: 3963 + ts - time stepping context 3964 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 3965 3966 Note: 3967 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 3968 3969 Level: intermediate 3970 3971 .seealso: TSGetCFLTime(), TSADAPTCFL 3972 @*/ 3973 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 3974 { 3975 PetscFunctionBegin; 3976 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3977 ts->cfltime_local = cfltime; 3978 ts->cfltime = -1.; 3979 PetscFunctionReturn(0); 3980 } 3981 3982 #undef __FUNCT__ 3983 #define __FUNCT__ "TSGetCFLTime" 3984 /*@ 3985 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 3986 3987 Collective on TS 3988 3989 Input Arguments: 3990 . ts - time stepping context 3991 3992 Output Arguments: 3993 . cfltime - maximum stable time step for forward Euler 3994 3995 Level: advanced 3996 3997 .seealso: TSSetCFLTimeLocal() 3998 @*/ 3999 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 4000 { 4001 PetscErrorCode ierr; 4002 4003 PetscFunctionBegin; 4004 if (ts->cfltime < 0) { 4005 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,((PetscObject)ts)->comm);CHKERRQ(ierr); 4006 } 4007 *cfltime = ts->cfltime; 4008 PetscFunctionReturn(0); 4009 } 4010 4011 #undef __FUNCT__ 4012 #define __FUNCT__ "TSVISetVariableBounds" 4013 /*@ 4014 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 4015 4016 Input Parameters: 4017 . ts - the TS context. 4018 . xl - lower bound. 4019 . xu - upper bound. 4020 4021 Notes: 4022 If this routine is not called then the lower and upper bounds are set to 4023 SNES_VI_NINF and SNES_VI_INF respectively during SNESSetUp(). 4024 4025 Level: advanced 4026 4027 @*/ 4028 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 4029 { 4030 PetscErrorCode ierr; 4031 SNES snes; 4032 4033 PetscFunctionBegin; 4034 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4035 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 4036 PetscFunctionReturn(0); 4037 } 4038 4039 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4040 #include <mex.h> 4041 4042 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 4043 4044 #undef __FUNCT__ 4045 #define __FUNCT__ "TSComputeFunction_Matlab" 4046 /* 4047 TSComputeFunction_Matlab - Calls the function that has been set with 4048 TSSetFunctionMatlab(). 4049 4050 Collective on TS 4051 4052 Input Parameters: 4053 + snes - the TS context 4054 - u - input vector 4055 4056 Output Parameter: 4057 . y - function vector, as set by TSSetFunction() 4058 4059 Notes: 4060 TSComputeFunction() is typically used within nonlinear solvers 4061 implementations, so most users would not generally call this routine 4062 themselves. 4063 4064 Level: developer 4065 4066 .keywords: TS, nonlinear, compute, function 4067 4068 .seealso: TSSetFunction(), TSGetFunction() 4069 */ 4070 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx) 4071 { 4072 PetscErrorCode ierr; 4073 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4074 int nlhs = 1,nrhs = 7; 4075 mxArray *plhs[1],*prhs[7]; 4076 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 4077 4078 PetscFunctionBegin; 4079 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 4080 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4081 PetscValidHeaderSpecific(udot,VEC_CLASSID,4); 4082 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 4083 PetscCheckSameComm(snes,1,u,3); 4084 PetscCheckSameComm(snes,1,y,5); 4085 4086 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4087 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4088 ierr = PetscMemcpy(&lxdot,&udot,sizeof(udot));CHKERRQ(ierr); 4089 ierr = PetscMemcpy(&ly,&y,sizeof(u));CHKERRQ(ierr); 4090 4091 prhs[0] = mxCreateDoubleScalar((double)ls); 4092 prhs[1] = mxCreateDoubleScalar(time); 4093 prhs[2] = mxCreateDoubleScalar((double)lx); 4094 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4095 prhs[4] = mxCreateDoubleScalar((double)ly); 4096 prhs[5] = mxCreateString(sctx->funcname); 4097 prhs[6] = sctx->ctx; 4098 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 4099 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4100 mxDestroyArray(prhs[0]); 4101 mxDestroyArray(prhs[1]); 4102 mxDestroyArray(prhs[2]); 4103 mxDestroyArray(prhs[3]); 4104 mxDestroyArray(prhs[4]); 4105 mxDestroyArray(prhs[5]); 4106 mxDestroyArray(plhs[0]); 4107 PetscFunctionReturn(0); 4108 } 4109 4110 4111 #undef __FUNCT__ 4112 #define __FUNCT__ "TSSetFunctionMatlab" 4113 /* 4114 TSSetFunctionMatlab - Sets the function evaluation routine and function 4115 vector for use by the TS routines in solving ODEs 4116 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4117 4118 Logically Collective on TS 4119 4120 Input Parameters: 4121 + ts - the TS context 4122 - func - function evaluation routine 4123 4124 Calling sequence of func: 4125 $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx); 4126 4127 Level: beginner 4128 4129 .keywords: TS, nonlinear, set, function 4130 4131 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4132 */ 4133 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 4134 { 4135 PetscErrorCode ierr; 4136 TSMatlabContext *sctx; 4137 4138 PetscFunctionBegin; 4139 /* currently sctx is memory bleed */ 4140 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4141 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4142 /* 4143 This should work, but it doesn't 4144 sctx->ctx = ctx; 4145 mexMakeArrayPersistent(sctx->ctx); 4146 */ 4147 sctx->ctx = mxDuplicateArray(ctx); 4148 4149 ierr = TSSetIFunction(ts,PETSC_NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4150 PetscFunctionReturn(0); 4151 } 4152 4153 #undef __FUNCT__ 4154 #define __FUNCT__ "TSComputeJacobian_Matlab" 4155 /* 4156 TSComputeJacobian_Matlab - Calls the function that has been set with 4157 TSSetJacobianMatlab(). 4158 4159 Collective on TS 4160 4161 Input Parameters: 4162 + ts - the TS context 4163 . u - input vector 4164 . A, B - the matrices 4165 - ctx - user context 4166 4167 Output Parameter: 4168 . flag - structure of the matrix 4169 4170 Level: developer 4171 4172 .keywords: TS, nonlinear, compute, function 4173 4174 .seealso: TSSetFunction(), TSGetFunction() 4175 @*/ 4176 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4177 { 4178 PetscErrorCode ierr; 4179 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4180 int nlhs = 2,nrhs = 9; 4181 mxArray *plhs[2],*prhs[9]; 4182 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 4183 4184 PetscFunctionBegin; 4185 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4186 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4187 4188 /* call Matlab function in ctx with arguments u and y */ 4189 4190 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4191 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4192 ierr = PetscMemcpy(&lxdot,&udot,sizeof(u));CHKERRQ(ierr); 4193 ierr = PetscMemcpy(&lA,A,sizeof(u));CHKERRQ(ierr); 4194 ierr = PetscMemcpy(&lB,B,sizeof(u));CHKERRQ(ierr); 4195 4196 prhs[0] = mxCreateDoubleScalar((double)ls); 4197 prhs[1] = mxCreateDoubleScalar((double)time); 4198 prhs[2] = mxCreateDoubleScalar((double)lx); 4199 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4200 prhs[4] = mxCreateDoubleScalar((double)shift); 4201 prhs[5] = mxCreateDoubleScalar((double)lA); 4202 prhs[6] = mxCreateDoubleScalar((double)lB); 4203 prhs[7] = mxCreateString(sctx->funcname); 4204 prhs[8] = sctx->ctx; 4205 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 4206 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4207 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4208 mxDestroyArray(prhs[0]); 4209 mxDestroyArray(prhs[1]); 4210 mxDestroyArray(prhs[2]); 4211 mxDestroyArray(prhs[3]); 4212 mxDestroyArray(prhs[4]); 4213 mxDestroyArray(prhs[5]); 4214 mxDestroyArray(prhs[6]); 4215 mxDestroyArray(prhs[7]); 4216 mxDestroyArray(plhs[0]); 4217 mxDestroyArray(plhs[1]); 4218 PetscFunctionReturn(0); 4219 } 4220 4221 4222 #undef __FUNCT__ 4223 #define __FUNCT__ "TSSetJacobianMatlab" 4224 /* 4225 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4226 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 4227 4228 Logically Collective on TS 4229 4230 Input Parameters: 4231 + ts - the TS context 4232 . A,B - Jacobian matrices 4233 . func - function evaluation routine 4234 - ctx - user context 4235 4236 Calling sequence of func: 4237 $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx); 4238 4239 4240 Level: developer 4241 4242 .keywords: TS, nonlinear, set, function 4243 4244 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4245 */ 4246 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 4247 { 4248 PetscErrorCode ierr; 4249 TSMatlabContext *sctx; 4250 4251 PetscFunctionBegin; 4252 /* currently sctx is memory bleed */ 4253 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4254 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4255 /* 4256 This should work, but it doesn't 4257 sctx->ctx = ctx; 4258 mexMakeArrayPersistent(sctx->ctx); 4259 */ 4260 sctx->ctx = mxDuplicateArray(ctx); 4261 4262 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4263 PetscFunctionReturn(0); 4264 } 4265 4266 #undef __FUNCT__ 4267 #define __FUNCT__ "TSMonitor_Matlab" 4268 /* 4269 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 4270 4271 Collective on TS 4272 4273 .seealso: TSSetFunction(), TSGetFunction() 4274 @*/ 4275 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx) 4276 { 4277 PetscErrorCode ierr; 4278 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4279 int nlhs = 1,nrhs = 6; 4280 mxArray *plhs[1],*prhs[6]; 4281 long long int lx = 0,ls = 0; 4282 4283 PetscFunctionBegin; 4284 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4285 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 4286 4287 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4288 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4289 4290 prhs[0] = mxCreateDoubleScalar((double)ls); 4291 prhs[1] = mxCreateDoubleScalar((double)it); 4292 prhs[2] = mxCreateDoubleScalar((double)time); 4293 prhs[3] = mxCreateDoubleScalar((double)lx); 4294 prhs[4] = mxCreateString(sctx->funcname); 4295 prhs[5] = sctx->ctx; 4296 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 4297 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4298 mxDestroyArray(prhs[0]); 4299 mxDestroyArray(prhs[1]); 4300 mxDestroyArray(prhs[2]); 4301 mxDestroyArray(prhs[3]); 4302 mxDestroyArray(prhs[4]); 4303 mxDestroyArray(plhs[0]); 4304 PetscFunctionReturn(0); 4305 } 4306 4307 4308 #undef __FUNCT__ 4309 #define __FUNCT__ "TSMonitorSetMatlab" 4310 /* 4311 TSMonitorSetMatlab - Sets the monitor function from Matlab 4312 4313 Level: developer 4314 4315 .keywords: TS, nonlinear, set, function 4316 4317 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4318 */ 4319 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 4320 { 4321 PetscErrorCode ierr; 4322 TSMatlabContext *sctx; 4323 4324 PetscFunctionBegin; 4325 /* currently sctx is memory bleed */ 4326 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4327 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4328 /* 4329 This should work, but it doesn't 4330 sctx->ctx = ctx; 4331 mexMakeArrayPersistent(sctx->ctx); 4332 */ 4333 sctx->ctx = mxDuplicateArray(ctx); 4334 4335 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 4336 PetscFunctionReturn(0); 4337 } 4338 #endif 4339 4340 4341 4342 #undef __FUNCT__ 4343 #define __FUNCT__ "TSMonitorLGSolution" 4344 /*@C 4345 TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector 4346 in a time based line graph 4347 4348 Collective on TS 4349 4350 Input Parameters: 4351 + ts - the TS context 4352 . step - current time-step 4353 . ptime - current time 4354 - lg - a line graph object 4355 4356 Level: intermediate 4357 4358 Notes: each process in a parallel run displays its component solutions in a separate window 4359 4360 .keywords: TS, vector, monitor, view 4361 4362 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4363 @*/ 4364 PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4365 { 4366 PetscErrorCode ierr; 4367 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4368 const PetscScalar *yy; 4369 PetscInt dim; 4370 4371 PetscFunctionBegin; 4372 if (!step) { 4373 PetscDrawAxis axis; 4374 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4375 ierr = PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");CHKERRQ(ierr); 4376 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4377 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4378 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4379 } 4380 ierr = VecGetArrayRead(u,&yy);CHKERRQ(ierr); 4381 #if defined(PETSC_USE_COMPLEX) 4382 { 4383 PetscReal *yreal; 4384 PetscInt i,n; 4385 ierr = VecGetLocalSize(u,&n);CHKERRQ(ierr); 4386 ierr = PetscMalloc(n*sizeof(PetscReal),&yreal);CHKERRQ(ierr); 4387 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4388 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4389 ierr = PetscFree(yreal);CHKERRQ(ierr); 4390 } 4391 #else 4392 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4393 #endif 4394 ierr = VecRestoreArrayRead(u,&yy);CHKERRQ(ierr); 4395 if (((ctx->howoften > 0) && (!(step % ctx->howoften)) && (step > -1)) || ((ctx->howoften == -1) && (step == -1))) { 4396 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4397 } 4398 PetscFunctionReturn(0); 4399 } 4400 4401 #undef __FUNCT__ 4402 #define __FUNCT__ "TSMonitorLGError" 4403 /*@C 4404 TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector 4405 in a time based line graph 4406 4407 Collective on TS 4408 4409 Input Parameters: 4410 + ts - the TS context 4411 . step - current time-step 4412 . ptime - current time 4413 - lg - a line graph object 4414 4415 Level: intermediate 4416 4417 Notes: 4418 Only for sequential solves. 4419 4420 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 4421 4422 Options Database Keys: 4423 . -ts_monitor_lg_error - create a graphical monitor of error history 4424 4425 .keywords: TS, vector, monitor, view 4426 4427 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 4428 @*/ 4429 PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4430 { 4431 PetscErrorCode ierr; 4432 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4433 const PetscScalar *yy; 4434 Vec y; 4435 PetscInt dim; 4436 4437 PetscFunctionBegin; 4438 if (!step) { 4439 PetscDrawAxis axis; 4440 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4441 ierr = PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");CHKERRQ(ierr); 4442 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4443 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4444 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4445 } 4446 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 4447 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 4448 ierr = VecAXPY(y,-1.0,u);CHKERRQ(ierr); 4449 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 4450 #if defined(PETSC_USE_COMPLEX) 4451 { 4452 PetscReal *yreal; 4453 PetscInt i,n; 4454 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 4455 ierr = PetscMalloc(n*sizeof(PetscReal),&yreal);CHKERRQ(ierr); 4456 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4457 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4458 ierr = PetscFree(yreal);CHKERRQ(ierr); 4459 } 4460 #else 4461 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4462 #endif 4463 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 4464 ierr = VecDestroy(&y);CHKERRQ(ierr); 4465 if (((ctx->howoften > 0) && (!(step % ctx->howoften)) && (step > -1)) || ((ctx->howoften == -1) && (step == -1))) { 4466 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4467 } 4468 PetscFunctionReturn(0); 4469 } 4470 4471 #undef __FUNCT__ 4472 #define __FUNCT__ "TSMonitorLGSNESIterations" 4473 PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4474 { 4475 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4476 PetscReal x = ptime,y; 4477 PetscErrorCode ierr; 4478 PetscInt its; 4479 4480 PetscFunctionBegin; 4481 if (!n) { 4482 PetscDrawAxis axis; 4483 4484 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4485 ierr = PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");CHKERRQ(ierr); 4486 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4487 4488 ctx->snes_its = 0; 4489 } 4490 ierr = TSGetSNESIterations(ts,&its);CHKERRQ(ierr); 4491 y = its - ctx->snes_its; 4492 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4493 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 4494 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4495 } 4496 ctx->snes_its = its; 4497 PetscFunctionReturn(0); 4498 } 4499 4500 #undef __FUNCT__ 4501 #define __FUNCT__ "TSMonitorLGKSPIterations" 4502 PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4503 { 4504 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4505 PetscReal x = ptime,y; 4506 PetscErrorCode ierr; 4507 PetscInt its; 4508 4509 PetscFunctionBegin; 4510 if (!n) { 4511 PetscDrawAxis axis; 4512 4513 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4514 ierr = PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");CHKERRQ(ierr); 4515 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4516 4517 ctx->ksp_its = 0; 4518 } 4519 ierr = TSGetKSPIterations(ts,&its);CHKERRQ(ierr); 4520 y = its - ctx->ksp_its; 4521 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4522 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 4523 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4524 } 4525 ctx->ksp_its = its; 4526 PetscFunctionReturn(0); 4527 } 4528 4529 #undef __FUNCT__ 4530 #define __FUNCT__ "TSComputeLinearStability" 4531 /*@ 4532 TSComputeLinearStability - computes the linear stability function at a point 4533 4534 Collective on TS and Vec 4535 4536 Input Parameters: 4537 + ts - the TS context 4538 - xr,xi - real and imaginary part of input arguments 4539 4540 Output Parameters: 4541 . yr,yi - real and imaginary part of function value 4542 4543 Level: developer 4544 4545 .keywords: TS, compute 4546 4547 .seealso: TSSetRHSFunction(), TSComputeIFunction() 4548 @*/ 4549 PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 4550 { 4551 PetscErrorCode ierr; 4552 4553 PetscFunctionBegin; 4554 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4555 if (!ts->ops->linearstability) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_SUP,"Linearized stability function not provided for this method"); 4556 ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr); 4557 PetscFunctionReturn(0); 4558 } 4559