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