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