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