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 (rhsjacobianfunc) { 364 ierr = PetscLogEventBegin(TS_JacobianEval,ts,U,*A,*B);CHKERRQ(ierr); 365 *flg = DIFFERENT_NONZERO_PATTERN; 366 PetscStackPush("TS user Jacobian function"); 367 ierr = (*rhsjacobianfunc)(ts,t,U,A,B,flg,ctx);CHKERRQ(ierr); 368 PetscStackPop; 369 ierr = PetscLogEventEnd(TS_JacobianEval,ts,U,*A,*B);CHKERRQ(ierr); 370 /* make sure user returned a correct Jacobian and preconditioner */ 371 PetscValidHeaderSpecific(*A,MAT_CLASSID,4); 372 PetscValidHeaderSpecific(*B,MAT_CLASSID,5); 373 } else { 374 ierr = MatZeroEntries(*A);CHKERRQ(ierr); 375 if (*A != *B) {ierr = MatZeroEntries(*B);CHKERRQ(ierr);} 376 *flg = SAME_NONZERO_PATTERN; 377 } 378 ts->rhsjacobian.time = t; 379 ts->rhsjacobian.X = U; 380 ierr = PetscObjectStateQuery((PetscObject)U,&ts->rhsjacobian.Xstate);CHKERRQ(ierr); 381 ts->rhsjacobian.mstructure = *flg; 382 PetscFunctionReturn(0); 383 } 384 385 #undef __FUNCT__ 386 #define __FUNCT__ "TSComputeRHSFunction" 387 /*@ 388 TSComputeRHSFunction - Evaluates the right-hand-side function. 389 390 Collective on TS and Vec 391 392 Input Parameters: 393 + ts - the TS context 394 . t - current time 395 - U - state vector 396 397 Output Parameter: 398 . y - right hand side 399 400 Note: 401 Most users should not need to explicitly call this routine, as it 402 is used internally within the nonlinear solvers. 403 404 Level: developer 405 406 .keywords: TS, compute 407 408 .seealso: TSSetRHSFunction(), TSComputeIFunction() 409 @*/ 410 PetscErrorCode TSComputeRHSFunction(TS ts,PetscReal t,Vec U,Vec y) 411 { 412 PetscErrorCode ierr; 413 TSRHSFunction rhsfunction; 414 TSIFunction ifunction; 415 void *ctx; 416 DM dm; 417 418 PetscFunctionBegin; 419 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 420 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 421 PetscValidHeaderSpecific(y,VEC_CLASSID,4); 422 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 423 ierr = DMTSGetRHSFunction(dm,&rhsfunction,&ctx);CHKERRQ(ierr); 424 ierr = DMTSGetIFunction(dm,&ifunction,NULL);CHKERRQ(ierr); 425 426 if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()"); 427 428 ierr = PetscLogEventBegin(TS_FunctionEval,ts,U,y,0);CHKERRQ(ierr); 429 if (rhsfunction) { 430 PetscStackPush("TS user right-hand-side function"); 431 ierr = (*rhsfunction)(ts,t,U,y,ctx);CHKERRQ(ierr); 432 PetscStackPop; 433 } else { 434 ierr = VecZeroEntries(y);CHKERRQ(ierr); 435 } 436 437 ierr = PetscLogEventEnd(TS_FunctionEval,ts,U,y,0);CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 } 440 441 #undef __FUNCT__ 442 #define __FUNCT__ "TSComputeSolutionFunction" 443 /*@ 444 TSComputeSolutionFunction - Evaluates the solution function. 445 446 Collective on TS and Vec 447 448 Input Parameters: 449 + ts - the TS context 450 - t - current time 451 452 Output Parameter: 453 . U - the solution 454 455 Note: 456 Most users should not need to explicitly call this routine, as it 457 is used internally within the nonlinear solvers. 458 459 Level: developer 460 461 .keywords: TS, compute 462 463 .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction() 464 @*/ 465 PetscErrorCode TSComputeSolutionFunction(TS ts,PetscReal t,Vec U) 466 { 467 PetscErrorCode ierr; 468 TSSolutionFunction solutionfunction; 469 void *ctx; 470 DM dm; 471 472 PetscFunctionBegin; 473 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 474 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 475 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 476 ierr = DMTSGetSolutionFunction(dm,&solutionfunction,&ctx);CHKERRQ(ierr); 477 478 if (solutionfunction) { 479 PetscStackPush("TS user solution function"); 480 ierr = (*solutionfunction)(ts,t,U,ctx);CHKERRQ(ierr); 481 PetscStackPop; 482 } 483 PetscFunctionReturn(0); 484 } 485 #undef __FUNCT__ 486 #define __FUNCT__ "TSComputeForcingFunction" 487 /*@ 488 TSComputeForcingFunction - Evaluates the forcing function. 489 490 Collective on TS and Vec 491 492 Input Parameters: 493 + ts - the TS context 494 - t - current time 495 496 Output Parameter: 497 . U - the function value 498 499 Note: 500 Most users should not need to explicitly call this routine, as it 501 is used internally within the nonlinear solvers. 502 503 Level: developer 504 505 .keywords: TS, compute 506 507 .seealso: TSSetSolutionFunction(), TSSetRHSFunction(), TSComputeIFunction() 508 @*/ 509 PetscErrorCode TSComputeForcingFunction(TS ts,PetscReal t,Vec U) 510 { 511 PetscErrorCode ierr, (*forcing)(TS,PetscReal,Vec,void*); 512 void *ctx; 513 DM dm; 514 515 PetscFunctionBegin; 516 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 517 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 518 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 519 ierr = DMTSGetForcingFunction(dm,&forcing,&ctx);CHKERRQ(ierr); 520 521 if (forcing) { 522 PetscStackPush("TS user forcing function"); 523 ierr = (*forcing)(ts,t,U,ctx);CHKERRQ(ierr); 524 PetscStackPop; 525 } 526 PetscFunctionReturn(0); 527 } 528 529 #undef __FUNCT__ 530 #define __FUNCT__ "TSGetRHSVec_Private" 531 static PetscErrorCode TSGetRHSVec_Private(TS ts,Vec *Frhs) 532 { 533 Vec F; 534 PetscErrorCode ierr; 535 536 PetscFunctionBegin; 537 *Frhs = NULL; 538 ierr = TSGetIFunction(ts,&F,NULL,NULL);CHKERRQ(ierr); 539 if (!ts->Frhs) { 540 ierr = VecDuplicate(F,&ts->Frhs);CHKERRQ(ierr); 541 } 542 *Frhs = ts->Frhs; 543 PetscFunctionReturn(0); 544 } 545 546 #undef __FUNCT__ 547 #define __FUNCT__ "TSGetRHSMats_Private" 548 static PetscErrorCode TSGetRHSMats_Private(TS ts,Mat *Arhs,Mat *Brhs) 549 { 550 Mat A,B; 551 PetscErrorCode ierr; 552 553 PetscFunctionBegin; 554 ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr); 555 if (Arhs) { 556 if (!ts->Arhs) { 557 ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&ts->Arhs);CHKERRQ(ierr); 558 } 559 *Arhs = ts->Arhs; 560 } 561 if (Brhs) { 562 if (!ts->Brhs) { 563 ierr = MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ts->Brhs);CHKERRQ(ierr); 564 } 565 *Brhs = ts->Brhs; 566 } 567 PetscFunctionReturn(0); 568 } 569 570 #undef __FUNCT__ 571 #define __FUNCT__ "TSComputeIFunction" 572 /*@ 573 TSComputeIFunction - Evaluates the DAE residual written in implicit form F(t,U,Udot)=0 574 575 Collective on TS and Vec 576 577 Input Parameters: 578 + ts - the TS context 579 . t - current time 580 . U - state vector 581 . Udot - time derivative of state vector 582 - imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate 583 584 Output Parameter: 585 . Y - right hand side 586 587 Note: 588 Most users should not need to explicitly call this routine, as it 589 is used internally within the nonlinear solvers. 590 591 If the user did did not write their equations in implicit form, this 592 function recasts them in implicit form. 593 594 Level: developer 595 596 .keywords: TS, compute 597 598 .seealso: TSSetIFunction(), TSComputeRHSFunction() 599 @*/ 600 PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec Y,PetscBool imex) 601 { 602 PetscErrorCode ierr; 603 TSIFunction ifunction; 604 TSRHSFunction rhsfunction; 605 void *ctx; 606 DM dm; 607 608 PetscFunctionBegin; 609 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 610 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 611 PetscValidHeaderSpecific(Udot,VEC_CLASSID,4); 612 PetscValidHeaderSpecific(Y,VEC_CLASSID,5); 613 614 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 615 ierr = DMTSGetIFunction(dm,&ifunction,&ctx);CHKERRQ(ierr); 616 ierr = DMTSGetRHSFunction(dm,&rhsfunction,NULL);CHKERRQ(ierr); 617 618 if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()"); 619 620 ierr = PetscLogEventBegin(TS_FunctionEval,ts,U,Udot,Y);CHKERRQ(ierr); 621 if (ifunction) { 622 PetscStackPush("TS user implicit function"); 623 ierr = (*ifunction)(ts,t,U,Udot,Y,ctx);CHKERRQ(ierr); 624 PetscStackPop; 625 } 626 if (imex) { 627 if (!ifunction) { 628 ierr = VecCopy(Udot,Y);CHKERRQ(ierr); 629 } 630 } else if (rhsfunction) { 631 if (ifunction) { 632 Vec Frhs; 633 ierr = TSGetRHSVec_Private(ts,&Frhs);CHKERRQ(ierr); 634 ierr = TSComputeRHSFunction(ts,t,U,Frhs);CHKERRQ(ierr); 635 ierr = VecAXPY(Y,-1,Frhs);CHKERRQ(ierr); 636 } else { 637 ierr = TSComputeRHSFunction(ts,t,U,Y);CHKERRQ(ierr); 638 ierr = VecAYPX(Y,-1,Udot);CHKERRQ(ierr); 639 } 640 } 641 ierr = PetscLogEventEnd(TS_FunctionEval,ts,U,Udot,Y);CHKERRQ(ierr); 642 PetscFunctionReturn(0); 643 } 644 645 #undef __FUNCT__ 646 #define __FUNCT__ "TSComputeIJacobian" 647 /*@ 648 TSComputeIJacobian - Evaluates the Jacobian of the DAE 649 650 Collective on TS and Vec 651 652 Input 653 Input Parameters: 654 + ts - the TS context 655 . t - current timestep 656 . U - state vector 657 . Udot - time derivative of state vector 658 . shift - shift to apply, see note below 659 - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate 660 661 Output Parameters: 662 + A - Jacobian matrix 663 . B - optional preconditioning matrix 664 - flag - flag indicating matrix structure 665 666 Notes: 667 If F(t,U,Udot)=0 is the DAE, the required Jacobian is 668 669 dF/dU + shift*dF/dUdot 670 671 Most users should not need to explicitly call this routine, as it 672 is used internally within the nonlinear solvers. 673 674 Level: developer 675 676 .keywords: TS, compute, Jacobian, matrix 677 678 .seealso: TSSetIJacobian() 679 @*/ 680 PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,PetscBool imex) 681 { 682 PetscInt Ustate, Udotstate; 683 PetscErrorCode ierr; 684 TSIJacobian ijacobian; 685 TSRHSJacobian rhsjacobian; 686 DM dm; 687 void *ctx; 688 689 PetscFunctionBegin; 690 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 691 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 692 PetscValidHeaderSpecific(Udot,VEC_CLASSID,4); 693 PetscValidPointer(A,6); 694 PetscValidHeaderSpecific(*A,MAT_CLASSID,6); 695 PetscValidPointer(B,7); 696 PetscValidHeaderSpecific(*B,MAT_CLASSID,7); 697 PetscValidPointer(flg,8); 698 699 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 700 ierr = DMTSGetIJacobian(dm,&ijacobian,&ctx);CHKERRQ(ierr); 701 ierr = DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);CHKERRQ(ierr); 702 703 ierr = PetscObjectStateQuery((PetscObject)U,&Ustate);CHKERRQ(ierr); 704 ierr = PetscObjectStateQuery((PetscObject)Udot,&Udotstate);CHKERRQ(ierr); 705 if (ts->ijacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->ijacobian.X == U && ts->ijacobian.Xstate == Ustate && ts->ijacobian.Xdot == Udot && ts->ijacobian.Xdotstate == Udotstate && ts->ijacobian.imex == imex))) { 706 *flg = ts->ijacobian.mstructure; 707 ierr = MatScale(*A, shift / ts->ijacobian.shift);CHKERRQ(ierr); 708 PetscFunctionReturn(0); 709 } 710 711 if (!rhsjacobian && !ijacobian) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()"); 712 713 *flg = SAME_NONZERO_PATTERN; /* In case we're solving a linear problem in which case it wouldn't get initialized below. */ 714 ierr = PetscLogEventBegin(TS_JacobianEval,ts,U,*A,*B);CHKERRQ(ierr); 715 if (ijacobian) { 716 *flg = DIFFERENT_NONZERO_PATTERN; 717 PetscStackPush("TS user implicit Jacobian"); 718 ierr = (*ijacobian)(ts,t,U,Udot,shift,A,B,flg,ctx);CHKERRQ(ierr); 719 PetscStackPop; 720 /* make sure user returned a correct Jacobian and preconditioner */ 721 PetscValidHeaderSpecific(*A,MAT_CLASSID,4); 722 PetscValidHeaderSpecific(*B,MAT_CLASSID,5); 723 } 724 if (imex) { 725 if (!ijacobian) { /* system was written as Udot = G(t,U) */ 726 ierr = MatZeroEntries(*A);CHKERRQ(ierr); 727 ierr = MatShift(*A,shift);CHKERRQ(ierr); 728 if (*A != *B) { 729 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 730 ierr = MatShift(*B,shift);CHKERRQ(ierr); 731 } 732 *flg = SAME_PRECONDITIONER; 733 } 734 } else { 735 if (!ijacobian) { 736 ierr = TSComputeRHSJacobian(ts,t,U,A,B,flg);CHKERRQ(ierr); 737 ierr = MatScale(*A,-1);CHKERRQ(ierr); 738 ierr = MatShift(*A,shift);CHKERRQ(ierr); 739 if (*A != *B) { 740 ierr = MatScale(*B,-1);CHKERRQ(ierr); 741 ierr = MatShift(*B,shift);CHKERRQ(ierr); 742 } 743 } else if (rhsjacobian) { 744 Mat Arhs,Brhs; 745 MatStructure axpy,flg2 = DIFFERENT_NONZERO_PATTERN; 746 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 747 ierr = TSComputeRHSJacobian(ts,t,U,&Arhs,&Brhs,&flg2);CHKERRQ(ierr); 748 axpy = (*flg == flg2) ? SAME_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN; 749 ierr = MatAXPY(*A,-1,Arhs,axpy);CHKERRQ(ierr); 750 if (*A != *B) { 751 ierr = MatAXPY(*B,-1,Brhs,axpy);CHKERRQ(ierr); 752 } 753 *flg = PetscMin(*flg,flg2); 754 } 755 } 756 757 ts->ijacobian.time = t; 758 ts->ijacobian.X = U; 759 ts->ijacobian.Xdot = Udot; 760 761 ierr = PetscObjectStateQuery((PetscObject)U,&ts->ijacobian.Xstate);CHKERRQ(ierr); 762 ierr = PetscObjectStateQuery((PetscObject)Udot,&ts->ijacobian.Xdotstate);CHKERRQ(ierr); 763 764 ts->ijacobian.shift = shift; 765 ts->ijacobian.imex = imex; 766 ts->ijacobian.mstructure = *flg; 767 768 ierr = PetscLogEventEnd(TS_JacobianEval,ts,U,*A,*B);CHKERRQ(ierr); 769 PetscFunctionReturn(0); 770 } 771 772 #undef __FUNCT__ 773 #define __FUNCT__ "TSSetRHSFunction" 774 /*@C 775 TSSetRHSFunction - Sets the routine for evaluating the function, 776 where U_t = G(t,u). 777 778 Logically Collective on TS 779 780 Input Parameters: 781 + ts - the TS context obtained from TSCreate() 782 . r - vector to put the computed right hand side (or NULL to have it created) 783 . f - routine for evaluating the right-hand-side function 784 - ctx - [optional] user-defined context for private data for the 785 function evaluation routine (may be NULL) 786 787 Calling sequence of func: 788 $ func (TS ts,PetscReal t,Vec u,Vec F,void *ctx); 789 790 + t - current timestep 791 . u - input vector 792 . F - function vector 793 - ctx - [optional] user-defined function context 794 795 Level: beginner 796 797 .keywords: TS, timestep, set, right-hand-side, function 798 799 .seealso: TSSetRHSJacobian(), TSSetIJacobian() 800 @*/ 801 PetscErrorCode TSSetRHSFunction(TS ts,Vec r,PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),void *ctx) 802 { 803 PetscErrorCode ierr; 804 SNES snes; 805 Vec ralloc = NULL; 806 DM dm; 807 808 PetscFunctionBegin; 809 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 810 if (r) PetscValidHeaderSpecific(r,VEC_CLASSID,2); 811 812 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 813 ierr = DMTSSetRHSFunction(dm,f,ctx);CHKERRQ(ierr); 814 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 815 if (!r && !ts->dm && ts->vec_sol) { 816 ierr = VecDuplicate(ts->vec_sol,&ralloc);CHKERRQ(ierr); 817 r = ralloc; 818 } 819 ierr = SNESSetFunction(snes,r,SNESTSFormFunction,ts);CHKERRQ(ierr); 820 ierr = VecDestroy(&ralloc);CHKERRQ(ierr); 821 PetscFunctionReturn(0); 822 } 823 824 #undef __FUNCT__ 825 #define __FUNCT__ "TSSetSolutionFunction" 826 /*@C 827 TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE 828 829 Logically Collective on TS 830 831 Input Parameters: 832 + ts - the TS context obtained from TSCreate() 833 . f - routine for evaluating the solution 834 - ctx - [optional] user-defined context for private data for the 835 function evaluation routine (may be NULL) 836 837 Calling sequence of func: 838 $ func (TS ts,PetscReal t,Vec u,void *ctx); 839 840 + t - current timestep 841 . u - output vector 842 - ctx - [optional] user-defined function context 843 844 Notes: 845 This routine is used for testing accuracy of time integration schemes when you already know the solution. 846 If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to 847 create closed-form solutions with non-physical forcing terms. 848 849 For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history. 850 851 Level: beginner 852 853 .keywords: TS, timestep, set, right-hand-side, function 854 855 .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetForcingFunction() 856 @*/ 857 PetscErrorCode TSSetSolutionFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx) 858 { 859 PetscErrorCode ierr; 860 DM dm; 861 862 PetscFunctionBegin; 863 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 864 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 865 ierr = DMTSSetSolutionFunction(dm,f,ctx);CHKERRQ(ierr); 866 PetscFunctionReturn(0); 867 } 868 869 #undef __FUNCT__ 870 #define __FUNCT__ "TSSetForcingFunction" 871 /*@C 872 TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE 873 874 Logically Collective on TS 875 876 Input Parameters: 877 + ts - the TS context obtained from TSCreate() 878 . f - routine for evaluating the forcing function 879 - ctx - [optional] user-defined context for private data for the 880 function evaluation routine (may be NULL) 881 882 Calling sequence of func: 883 $ func (TS ts,PetscReal t,Vec u,void *ctx); 884 885 + t - current timestep 886 . u - output vector 887 - ctx - [optional] user-defined function context 888 889 Notes: 890 This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to 891 create closed-form solutions with a non-physical forcing term. 892 893 For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history. 894 895 Level: beginner 896 897 .keywords: TS, timestep, set, right-hand-side, function 898 899 .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetSolutionFunction() 900 @*/ 901 PetscErrorCode TSSetForcingFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx) 902 { 903 PetscErrorCode ierr; 904 DM dm; 905 906 PetscFunctionBegin; 907 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 908 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 909 ierr = DMTSSetForcingFunction(dm,f,ctx);CHKERRQ(ierr); 910 PetscFunctionReturn(0); 911 } 912 913 #undef __FUNCT__ 914 #define __FUNCT__ "TSSetRHSJacobian" 915 /*@C 916 TSSetRHSJacobian - Sets the function to compute the Jacobian of F, 917 where U_t = G(U,t), as well as the location to store the matrix. 918 919 Logically Collective on TS 920 921 Input Parameters: 922 + ts - the TS context obtained from TSCreate() 923 . Amat - (approximate) Jacobian matrix 924 . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat) 925 . f - the Jacobian evaluation routine 926 - ctx - [optional] user-defined context for private data for the 927 Jacobian evaluation routine (may be NULL) 928 929 Calling sequence of func: 930 $ func (TS ts,PetscReal t,Vec u,Mat *A,Mat *B,MatStructure *flag,void *ctx); 931 932 + t - current timestep 933 . u - input vector 934 . Amat - (approximate) Jacobian matrix 935 . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat) 936 . flag - flag indicating information about the preconditioner matrix 937 structure (same as flag in KSPSetOperators()) 938 - ctx - [optional] user-defined context for matrix evaluation routine 939 940 Notes: 941 See KSPSetOperators() for important information about setting the flag 942 output parameter in the routine func(). Be sure to read this information! 943 944 The routine func() takes Mat * as the matrix arguments rather than Mat. 945 This allows the matrix evaluation routine to replace A and/or B with a 946 completely new matrix structure (not just different matrix elements) 947 when appropriate, for instance, if the nonzero structure is changing 948 throughout the global iterations. 949 950 Level: beginner 951 952 .keywords: TS, timestep, set, right-hand-side, Jacobian 953 954 .seealso: SNESComputeJacobianDefaultColor(), TSSetRHSFunction() 955 956 @*/ 957 PetscErrorCode TSSetRHSJacobian(TS ts,Mat Amat,Mat Pmat,TSRHSJacobian f,void *ctx) 958 { 959 PetscErrorCode ierr; 960 SNES snes; 961 DM dm; 962 TSIJacobian ijacobian; 963 964 PetscFunctionBegin; 965 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 966 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 967 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 968 if (Amat) PetscCheckSameComm(ts,1,Amat,2); 969 if (Pmat) PetscCheckSameComm(ts,1,Pmat,3); 970 971 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 972 ierr = DMTSSetRHSJacobian(dm,f,ctx);CHKERRQ(ierr); 973 ierr = DMTSGetIJacobian(dm,&ijacobian,NULL);CHKERRQ(ierr); 974 975 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 976 if (!ijacobian) { 977 ierr = SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);CHKERRQ(ierr); 978 } 979 if (Amat) { 980 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 981 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 982 983 ts->Arhs = Amat; 984 } 985 if (Pmat) { 986 ierr = PetscObjectReference((PetscObject)Pmat);CHKERRQ(ierr); 987 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 988 989 ts->Brhs = Pmat; 990 } 991 PetscFunctionReturn(0); 992 } 993 994 995 #undef __FUNCT__ 996 #define __FUNCT__ "TSSetIFunction" 997 /*@C 998 TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved. 999 1000 Logically Collective on TS 1001 1002 Input Parameters: 1003 + ts - the TS context obtained from TSCreate() 1004 . r - vector to hold the residual (or NULL to have it created internally) 1005 . f - the function evaluation routine 1006 - ctx - user-defined context for private data for the function evaluation routine (may be NULL) 1007 1008 Calling sequence of f: 1009 $ f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx); 1010 1011 + t - time at step/stage being solved 1012 . u - state vector 1013 . u_t - time derivative of state vector 1014 . F - function vector 1015 - ctx - [optional] user-defined context for matrix evaluation routine 1016 1017 Important: 1018 The user MUST call either this routine, TSSetRHSFunction(). This routine must be used when not solving an ODE, for example a DAE. 1019 1020 Level: beginner 1021 1022 .keywords: TS, timestep, set, DAE, Jacobian 1023 1024 .seealso: TSSetRHSJacobian(), TSSetRHSFunction(), TSSetIJacobian() 1025 @*/ 1026 PetscErrorCode TSSetIFunction(TS ts,Vec res,TSIFunction f,void *ctx) 1027 { 1028 PetscErrorCode ierr; 1029 SNES snes; 1030 Vec resalloc = NULL; 1031 DM dm; 1032 1033 PetscFunctionBegin; 1034 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1035 if (res) PetscValidHeaderSpecific(res,VEC_CLASSID,2); 1036 1037 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1038 ierr = DMTSSetIFunction(dm,f,ctx);CHKERRQ(ierr); 1039 1040 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1041 if (!res && !ts->dm && ts->vec_sol) { 1042 ierr = VecDuplicate(ts->vec_sol,&resalloc);CHKERRQ(ierr); 1043 res = resalloc; 1044 } 1045 ierr = SNESSetFunction(snes,res,SNESTSFormFunction,ts);CHKERRQ(ierr); 1046 ierr = VecDestroy(&resalloc);CHKERRQ(ierr); 1047 PetscFunctionReturn(0); 1048 } 1049 1050 #undef __FUNCT__ 1051 #define __FUNCT__ "TSGetIFunction" 1052 /*@C 1053 TSGetIFunction - Returns the vector where the implicit residual is stored and the function/contex to compute it. 1054 1055 Not Collective 1056 1057 Input Parameter: 1058 . ts - the TS context 1059 1060 Output Parameter: 1061 + r - vector to hold residual (or NULL) 1062 . func - the function to compute residual (or NULL) 1063 - ctx - the function context (or NULL) 1064 1065 Level: advanced 1066 1067 .keywords: TS, nonlinear, get, function 1068 1069 .seealso: TSSetIFunction(), SNESGetFunction() 1070 @*/ 1071 PetscErrorCode TSGetIFunction(TS ts,Vec *r,TSIFunction *func,void **ctx) 1072 { 1073 PetscErrorCode ierr; 1074 SNES snes; 1075 DM dm; 1076 1077 PetscFunctionBegin; 1078 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1079 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1080 ierr = SNESGetFunction(snes,r,NULL,NULL);CHKERRQ(ierr); 1081 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1082 ierr = DMTSGetIFunction(dm,func,ctx);CHKERRQ(ierr); 1083 PetscFunctionReturn(0); 1084 } 1085 1086 #undef __FUNCT__ 1087 #define __FUNCT__ "TSGetRHSFunction" 1088 /*@C 1089 TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it. 1090 1091 Not Collective 1092 1093 Input Parameter: 1094 . ts - the TS context 1095 1096 Output Parameter: 1097 + r - vector to hold computed right hand side (or NULL) 1098 . func - the function to compute right hand side (or NULL) 1099 - ctx - the function context (or NULL) 1100 1101 Level: advanced 1102 1103 .keywords: TS, nonlinear, get, function 1104 1105 .seealso: TSSetRhsfunction(), SNESGetFunction() 1106 @*/ 1107 PetscErrorCode TSGetRHSFunction(TS ts,Vec *r,TSRHSFunction *func,void **ctx) 1108 { 1109 PetscErrorCode ierr; 1110 SNES snes; 1111 DM dm; 1112 1113 PetscFunctionBegin; 1114 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1115 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1116 ierr = SNESGetFunction(snes,r,NULL,NULL);CHKERRQ(ierr); 1117 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1118 ierr = DMTSGetRHSFunction(dm,func,ctx);CHKERRQ(ierr); 1119 PetscFunctionReturn(0); 1120 } 1121 1122 #undef __FUNCT__ 1123 #define __FUNCT__ "TSSetIJacobian" 1124 /*@C 1125 TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function 1126 you provided with TSSetIFunction(). 1127 1128 Logically Collective on TS 1129 1130 Input Parameters: 1131 + ts - the TS context obtained from TSCreate() 1132 . Amat - (approximate) Jacobian matrix 1133 . Pmat - matrix used to compute preconditioner (usually the same as Amat) 1134 . f - the Jacobian evaluation routine 1135 - ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL) 1136 1137 Calling sequence of f: 1138 $ f(TS ts,PetscReal t,Vec U,Vec U_t,PetscReal a,Mat *Amat,Mat *Pmat,MatStructure *flag,void *ctx); 1139 1140 + t - time at step/stage being solved 1141 . U - state vector 1142 . U_t - time derivative of state vector 1143 . a - shift 1144 . Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t 1145 . Pmat - matrix used for constructing preconditioner, usually the same as Amat 1146 . flag - flag indicating information about the preconditioner matrix 1147 structure (same as flag in KSPSetOperators()) 1148 - ctx - [optional] user-defined context for matrix evaluation routine 1149 1150 Notes: 1151 The matrices Amat and Pmat are exactly the matrices that are used by SNES for the nonlinear solve. 1152 1153 The matrix dF/dU + a*dF/dU_t you provide turns out to be 1154 the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved. 1155 The time integrator internally approximates U_t by W+a*U where the positive "shift" 1156 a and vector W depend on the integration method, step size, and past states. For example with 1157 the backward Euler method a = 1/dt and W = -a*U(previous timestep) so 1158 W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt 1159 1160 Level: beginner 1161 1162 .keywords: TS, timestep, DAE, Jacobian 1163 1164 .seealso: TSSetIFunction(), TSSetRHSJacobian(), SNESComputeJacobianDefaultColor(), SNESComputeJacobianDefault() 1165 1166 @*/ 1167 PetscErrorCode TSSetIJacobian(TS ts,Mat Amat,Mat Pmat,TSIJacobian f,void *ctx) 1168 { 1169 PetscErrorCode ierr; 1170 SNES snes; 1171 DM dm; 1172 1173 PetscFunctionBegin; 1174 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1175 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 1176 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 1177 if (Amat) PetscCheckSameComm(ts,1,Amat,2); 1178 if (Pmat) PetscCheckSameComm(ts,1,Pmat,3); 1179 1180 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1181 ierr = DMTSSetIJacobian(dm,f,ctx);CHKERRQ(ierr); 1182 1183 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1184 ierr = SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1185 PetscFunctionReturn(0); 1186 } 1187 1188 #undef __FUNCT__ 1189 #define __FUNCT__ "TSLoad" 1190 /*@C 1191 TSLoad - Loads a KSP that has been stored in binary with KSPView(). 1192 1193 Collective on PetscViewer 1194 1195 Input Parameters: 1196 + newdm - the newly loaded TS, this needs to have been created with TSCreate() or 1197 some related function before a call to TSLoad(). 1198 - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() 1199 1200 Level: intermediate 1201 1202 Notes: 1203 The type is determined by the data in the file, any type set into the TS before this call is ignored. 1204 1205 Notes for advanced users: 1206 Most users should not need to know the details of the binary storage 1207 format, since TSLoad() and TSView() completely hide these details. 1208 But for anyone who's interested, the standard binary matrix storage 1209 format is 1210 .vb 1211 has not yet been determined 1212 .ve 1213 1214 .seealso: PetscViewerBinaryOpen(), TSView(), MatLoad(), VecLoad() 1215 @*/ 1216 PetscErrorCode TSLoad(TS ts, PetscViewer viewer) 1217 { 1218 PetscErrorCode ierr; 1219 PetscBool isbinary; 1220 PetscInt classid; 1221 char type[256]; 1222 DMTS sdm; 1223 DM dm; 1224 1225 PetscFunctionBegin; 1226 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1227 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1228 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1229 if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()"); 1230 1231 ierr = PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr); 1232 if (classid != TS_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Not TS next in file"); 1233 ierr = PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);CHKERRQ(ierr); 1234 ierr = TSSetType(ts, type);CHKERRQ(ierr); 1235 if (ts->ops->load) { 1236 ierr = (*ts->ops->load)(ts,viewer);CHKERRQ(ierr); 1237 } 1238 ierr = DMCreate(PetscObjectComm((PetscObject)ts),&dm);CHKERRQ(ierr); 1239 ierr = DMLoad(dm,viewer);CHKERRQ(ierr); 1240 ierr = TSSetDM(ts,dm);CHKERRQ(ierr); 1241 ierr = DMCreateGlobalVector(ts->dm,&ts->vec_sol);CHKERRQ(ierr); 1242 ierr = VecLoad(ts->vec_sol,viewer);CHKERRQ(ierr); 1243 ierr = DMGetDMTS(ts->dm,&sdm);CHKERRQ(ierr); 1244 ierr = DMTSLoad(sdm,viewer);CHKERRQ(ierr); 1245 PetscFunctionReturn(0); 1246 } 1247 1248 #include <petscdraw.h> 1249 #if defined(PETSC_HAVE_AMS) 1250 #include <petscviewerams.h> 1251 #endif 1252 #undef __FUNCT__ 1253 #define __FUNCT__ "TSView" 1254 /*@C 1255 TSView - Prints the TS data structure. 1256 1257 Collective on TS 1258 1259 Input Parameters: 1260 + ts - the TS context obtained from TSCreate() 1261 - viewer - visualization context 1262 1263 Options Database Key: 1264 . -ts_view - calls TSView() at end of TSStep() 1265 1266 Notes: 1267 The available visualization contexts include 1268 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 1269 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 1270 output where only the first processor opens 1271 the file. All other processors send their 1272 data to the first processor to print. 1273 1274 The user can open an alternative visualization context with 1275 PetscViewerASCIIOpen() - output to a specified file. 1276 1277 Level: beginner 1278 1279 .keywords: TS, timestep, view 1280 1281 .seealso: PetscViewerASCIIOpen() 1282 @*/ 1283 PetscErrorCode TSView(TS ts,PetscViewer viewer) 1284 { 1285 PetscErrorCode ierr; 1286 TSType type; 1287 PetscBool iascii,isstring,isundials,isbinary,isdraw; 1288 DMTS sdm; 1289 #if defined(PETSC_HAVE_AMS) 1290 PetscBool isams; 1291 #endif 1292 1293 PetscFunctionBegin; 1294 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1295 if (!viewer) { 1296 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts),&viewer);CHKERRQ(ierr); 1297 } 1298 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1299 PetscCheckSameComm(ts,1,viewer,2); 1300 1301 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1302 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 1303 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1304 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1305 #if defined(PETSC_HAVE_AMS) 1306 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERAMS,&isams);CHKERRQ(ierr); 1307 #endif 1308 if (iascii) { 1309 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer,"TS Object");CHKERRQ(ierr); 1310 ierr = PetscViewerASCIIPrintf(viewer," maximum steps=%D\n",ts->max_steps);CHKERRQ(ierr); 1311 ierr = PetscViewerASCIIPrintf(viewer," maximum time=%G\n",ts->max_time);CHKERRQ(ierr); 1312 if (ts->problem_type == TS_NONLINEAR) { 1313 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solver iterations=%D\n",ts->snes_its);CHKERRQ(ierr); 1314 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solve failures=%D\n",ts->num_snes_failures);CHKERRQ(ierr); 1315 } 1316 ierr = PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",ts->ksp_its);CHKERRQ(ierr); 1317 ierr = PetscViewerASCIIPrintf(viewer," total number of rejected steps=%D\n",ts->reject);CHKERRQ(ierr); 1318 ierr = DMGetDMTS(ts->dm,&sdm);CHKERRQ(ierr); 1319 ierr = DMTSView(sdm,viewer);CHKERRQ(ierr); 1320 if (ts->ops->view) { 1321 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1322 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1323 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1324 } 1325 } else if (isstring) { 1326 ierr = TSGetType(ts,&type);CHKERRQ(ierr); 1327 ierr = PetscViewerStringSPrintf(viewer," %-7.7s",type);CHKERRQ(ierr); 1328 } else if (isbinary) { 1329 PetscInt classid = TS_FILE_CLASSID; 1330 MPI_Comm comm; 1331 PetscMPIInt rank; 1332 char type[256]; 1333 1334 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 1335 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1336 if (!rank) { 1337 ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1338 ierr = PetscStrncpy(type,((PetscObject)ts)->type_name,256);CHKERRQ(ierr); 1339 ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr); 1340 } 1341 if (ts->ops->view) { 1342 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1343 } 1344 ierr = DMView(ts->dm,viewer);CHKERRQ(ierr); 1345 ierr = VecView(ts->vec_sol,viewer);CHKERRQ(ierr); 1346 ierr = DMGetDMTS(ts->dm,&sdm);CHKERRQ(ierr); 1347 ierr = DMTSView(sdm,viewer);CHKERRQ(ierr); 1348 } else if (isdraw) { 1349 PetscDraw draw; 1350 char str[36]; 1351 PetscReal x,y,bottom,h; 1352 1353 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1354 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 1355 ierr = PetscStrcpy(str,"TS: ");CHKERRQ(ierr); 1356 ierr = PetscStrcat(str,((PetscObject)ts)->type_name);CHKERRQ(ierr); 1357 ierr = PetscDrawBoxedString(draw,x,y,PETSC_DRAW_BLACK,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 1358 bottom = y - h; 1359 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr); 1360 if (ts->ops->view) { 1361 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1362 } 1363 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 1364 #if defined(PETSC_HAVE_AMS) 1365 } else if (isams) { 1366 if (((PetscObject)ts)->amsmem == -1) { 1367 ierr = PetscObjectViewAMS((PetscObject)ts,viewer);CHKERRQ(ierr); 1368 PetscStackCallAMS(AMS_Memory_take_access,(((PetscObject)ts)->amsmem)); 1369 PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)ts)->amsmem,"time step",&ts->steps,1,AMS_INT,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF)); 1370 PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)ts)->amsmem,"time",&ts->ptime,1,AMS_DOUBLE,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF)); 1371 PetscStackCallAMS(AMS_Memory_grant_access,(((PetscObject)ts)->amsmem)); 1372 } 1373 if (ts->ops->view) { 1374 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1375 } 1376 #endif 1377 } 1378 1379 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1380 ierr = PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);CHKERRQ(ierr); 1381 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1382 PetscFunctionReturn(0); 1383 } 1384 1385 1386 #undef __FUNCT__ 1387 #define __FUNCT__ "TSSetApplicationContext" 1388 /*@ 1389 TSSetApplicationContext - Sets an optional user-defined context for 1390 the timesteppers. 1391 1392 Logically Collective on TS 1393 1394 Input Parameters: 1395 + ts - the TS context obtained from TSCreate() 1396 - usrP - optional user context 1397 1398 Level: intermediate 1399 1400 .keywords: TS, timestep, set, application, context 1401 1402 .seealso: TSGetApplicationContext() 1403 @*/ 1404 PetscErrorCode TSSetApplicationContext(TS ts,void *usrP) 1405 { 1406 PetscFunctionBegin; 1407 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1408 ts->user = usrP; 1409 PetscFunctionReturn(0); 1410 } 1411 1412 #undef __FUNCT__ 1413 #define __FUNCT__ "TSGetApplicationContext" 1414 /*@ 1415 TSGetApplicationContext - Gets the user-defined context for the 1416 timestepper. 1417 1418 Not Collective 1419 1420 Input Parameter: 1421 . ts - the TS context obtained from TSCreate() 1422 1423 Output Parameter: 1424 . usrP - user context 1425 1426 Level: intermediate 1427 1428 .keywords: TS, timestep, get, application, context 1429 1430 .seealso: TSSetApplicationContext() 1431 @*/ 1432 PetscErrorCode TSGetApplicationContext(TS ts,void *usrP) 1433 { 1434 PetscFunctionBegin; 1435 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1436 *(void**)usrP = ts->user; 1437 PetscFunctionReturn(0); 1438 } 1439 1440 #undef __FUNCT__ 1441 #define __FUNCT__ "TSGetTimeStepNumber" 1442 /*@ 1443 TSGetTimeStepNumber - Gets the number of time steps completed. 1444 1445 Not Collective 1446 1447 Input Parameter: 1448 . ts - the TS context obtained from TSCreate() 1449 1450 Output Parameter: 1451 . iter - number of steps completed so far 1452 1453 Level: intermediate 1454 1455 .keywords: TS, timestep, get, iteration, number 1456 .seealso: TSGetTime(), TSGetTimeStep(), TSSetPreStep(), TSSetPreStage(), TSSetPostStage(), TSSetPostStep() 1457 @*/ 1458 PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt *iter) 1459 { 1460 PetscFunctionBegin; 1461 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1462 PetscValidIntPointer(iter,2); 1463 *iter = ts->steps; 1464 PetscFunctionReturn(0); 1465 } 1466 1467 #undef __FUNCT__ 1468 #define __FUNCT__ "TSSetInitialTimeStep" 1469 /*@ 1470 TSSetInitialTimeStep - Sets the initial timestep to be used, 1471 as well as the initial time. 1472 1473 Logically Collective on TS 1474 1475 Input Parameters: 1476 + ts - the TS context obtained from TSCreate() 1477 . initial_time - the initial time 1478 - time_step - the size of the timestep 1479 1480 Level: intermediate 1481 1482 .seealso: TSSetTimeStep(), TSGetTimeStep() 1483 1484 .keywords: TS, set, initial, timestep 1485 @*/ 1486 PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step) 1487 { 1488 PetscErrorCode ierr; 1489 1490 PetscFunctionBegin; 1491 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1492 ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr); 1493 ierr = TSSetTime(ts,initial_time);CHKERRQ(ierr); 1494 PetscFunctionReturn(0); 1495 } 1496 1497 #undef __FUNCT__ 1498 #define __FUNCT__ "TSSetTimeStep" 1499 /*@ 1500 TSSetTimeStep - Allows one to reset the timestep at any time, 1501 useful for simple pseudo-timestepping codes. 1502 1503 Logically Collective on TS 1504 1505 Input Parameters: 1506 + ts - the TS context obtained from TSCreate() 1507 - time_step - the size of the timestep 1508 1509 Level: intermediate 1510 1511 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1512 1513 .keywords: TS, set, timestep 1514 @*/ 1515 PetscErrorCode TSSetTimeStep(TS ts,PetscReal time_step) 1516 { 1517 PetscFunctionBegin; 1518 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1519 PetscValidLogicalCollectiveReal(ts,time_step,2); 1520 ts->time_step = time_step; 1521 ts->time_step_orig = time_step; 1522 PetscFunctionReturn(0); 1523 } 1524 1525 #undef __FUNCT__ 1526 #define __FUNCT__ "TSSetExactFinalTime" 1527 /*@ 1528 TSSetExactFinalTime - Determines whether to adapt the final time step to 1529 match the exact final time, interpolate solution to the exact final time, 1530 or just return at the final time TS computed. 1531 1532 Logically Collective on TS 1533 1534 Input Parameter: 1535 + ts - the time-step context 1536 - eftopt - exact final time option 1537 1538 Level: beginner 1539 1540 .seealso: TSExactFinalTimeOption 1541 @*/ 1542 PetscErrorCode TSSetExactFinalTime(TS ts,TSExactFinalTimeOption eftopt) 1543 { 1544 PetscFunctionBegin; 1545 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1546 PetscValidLogicalCollectiveEnum(ts,eftopt,2); 1547 ts->exact_final_time = eftopt; 1548 PetscFunctionReturn(0); 1549 } 1550 1551 #undef __FUNCT__ 1552 #define __FUNCT__ "TSGetTimeStep" 1553 /*@ 1554 TSGetTimeStep - Gets the current timestep size. 1555 1556 Not Collective 1557 1558 Input Parameter: 1559 . ts - the TS context obtained from TSCreate() 1560 1561 Output Parameter: 1562 . dt - the current timestep size 1563 1564 Level: intermediate 1565 1566 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1567 1568 .keywords: TS, get, timestep 1569 @*/ 1570 PetscErrorCode TSGetTimeStep(TS ts,PetscReal *dt) 1571 { 1572 PetscFunctionBegin; 1573 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1574 PetscValidRealPointer(dt,2); 1575 *dt = ts->time_step; 1576 PetscFunctionReturn(0); 1577 } 1578 1579 #undef __FUNCT__ 1580 #define __FUNCT__ "TSGetSolution" 1581 /*@ 1582 TSGetSolution - Returns the solution at the present timestep. It 1583 is valid to call this routine inside the function that you are evaluating 1584 in order to move to the new timestep. This vector not changed until 1585 the solution at the next timestep has been calculated. 1586 1587 Not Collective, but Vec returned is parallel if TS is parallel 1588 1589 Input Parameter: 1590 . ts - the TS context obtained from TSCreate() 1591 1592 Output Parameter: 1593 . v - the vector containing the solution 1594 1595 Level: intermediate 1596 1597 .seealso: TSGetTimeStep() 1598 1599 .keywords: TS, timestep, get, solution 1600 @*/ 1601 PetscErrorCode TSGetSolution(TS ts,Vec *v) 1602 { 1603 PetscFunctionBegin; 1604 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1605 PetscValidPointer(v,2); 1606 *v = ts->vec_sol; 1607 PetscFunctionReturn(0); 1608 } 1609 1610 /* ----- Routines to initialize and destroy a timestepper ---- */ 1611 #undef __FUNCT__ 1612 #define __FUNCT__ "TSSetProblemType" 1613 /*@ 1614 TSSetProblemType - Sets the type of problem to be solved. 1615 1616 Not collective 1617 1618 Input Parameters: 1619 + ts - The TS 1620 - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1621 .vb 1622 U_t - A U = 0 (linear) 1623 U_t - A(t) U = 0 (linear) 1624 F(t,U,U_t) = 0 (nonlinear) 1625 .ve 1626 1627 Level: beginner 1628 1629 .keywords: TS, problem type 1630 .seealso: TSSetUp(), TSProblemType, TS 1631 @*/ 1632 PetscErrorCode TSSetProblemType(TS ts, TSProblemType type) 1633 { 1634 PetscErrorCode ierr; 1635 1636 PetscFunctionBegin; 1637 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1638 ts->problem_type = type; 1639 if (type == TS_LINEAR) { 1640 SNES snes; 1641 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1642 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1643 } 1644 PetscFunctionReturn(0); 1645 } 1646 1647 #undef __FUNCT__ 1648 #define __FUNCT__ "TSGetProblemType" 1649 /*@C 1650 TSGetProblemType - Gets the type of problem to be solved. 1651 1652 Not collective 1653 1654 Input Parameter: 1655 . ts - The TS 1656 1657 Output Parameter: 1658 . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1659 .vb 1660 M U_t = A U 1661 M(t) U_t = A(t) U 1662 F(t,U,U_t) 1663 .ve 1664 1665 Level: beginner 1666 1667 .keywords: TS, problem type 1668 .seealso: TSSetUp(), TSProblemType, TS 1669 @*/ 1670 PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type) 1671 { 1672 PetscFunctionBegin; 1673 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1674 PetscValidIntPointer(type,2); 1675 *type = ts->problem_type; 1676 PetscFunctionReturn(0); 1677 } 1678 1679 #undef __FUNCT__ 1680 #define __FUNCT__ "TSSetUp" 1681 /*@ 1682 TSSetUp - Sets up the internal data structures for the later use 1683 of a timestepper. 1684 1685 Collective on TS 1686 1687 Input Parameter: 1688 . ts - the TS context obtained from TSCreate() 1689 1690 Notes: 1691 For basic use of the TS solvers the user need not explicitly call 1692 TSSetUp(), since these actions will automatically occur during 1693 the call to TSStep(). However, if one wishes to control this 1694 phase separately, TSSetUp() should be called after TSCreate() 1695 and optional routines of the form TSSetXXX(), but before TSStep(). 1696 1697 Level: advanced 1698 1699 .keywords: TS, timestep, setup 1700 1701 .seealso: TSCreate(), TSStep(), TSDestroy() 1702 @*/ 1703 PetscErrorCode TSSetUp(TS ts) 1704 { 1705 PetscErrorCode ierr; 1706 DM dm; 1707 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 1708 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 1709 TSIJacobian ijac; 1710 TSRHSJacobian rhsjac; 1711 1712 PetscFunctionBegin; 1713 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1714 if (ts->setupcalled) PetscFunctionReturn(0); 1715 1716 if (!((PetscObject)ts)->type_name) { 1717 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 1718 } 1719 1720 if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first"); 1721 1722 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 1723 1724 if (ts->ops->setup) { 1725 ierr = (*ts->ops->setup)(ts);CHKERRQ(ierr); 1726 } 1727 1728 /* in the case where we've set a DMTSFunction or what have you, we need the default SNESFunction 1729 to be set right but can't do it elsewhere due to the overreliance on ctx=ts. 1730 */ 1731 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1732 ierr = DMSNESGetFunction(dm,&func,NULL);CHKERRQ(ierr); 1733 if (!func) { 1734 ierr =DMSNESSetFunction(dm,SNESTSFormFunction,ts);CHKERRQ(ierr); 1735 } 1736 /* if the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it. 1737 Otherwise, the SNES will use coloring internally to form the Jacobian. 1738 */ 1739 ierr = DMSNESGetJacobian(dm,&jac,NULL);CHKERRQ(ierr); 1740 ierr = DMTSGetIJacobian(dm,&ijac,NULL);CHKERRQ(ierr); 1741 ierr = DMTSGetRHSJacobian(dm,&rhsjac,NULL);CHKERRQ(ierr); 1742 if (!jac && (ijac || rhsjac)) { 1743 ierr = DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1744 } 1745 ts->setupcalled = PETSC_TRUE; 1746 PetscFunctionReturn(0); 1747 } 1748 1749 #undef __FUNCT__ 1750 #define __FUNCT__ "TSReset" 1751 /*@ 1752 TSReset - Resets a TS context and removes any allocated Vecs and Mats. 1753 1754 Collective on TS 1755 1756 Input Parameter: 1757 . ts - the TS context obtained from TSCreate() 1758 1759 Level: beginner 1760 1761 .keywords: TS, timestep, reset 1762 1763 .seealso: TSCreate(), TSSetup(), TSDestroy() 1764 @*/ 1765 PetscErrorCode TSReset(TS ts) 1766 { 1767 PetscErrorCode ierr; 1768 1769 PetscFunctionBegin; 1770 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1771 if (ts->ops->reset) { 1772 ierr = (*ts->ops->reset)(ts);CHKERRQ(ierr); 1773 } 1774 if (ts->snes) {ierr = SNESReset(ts->snes);CHKERRQ(ierr);} 1775 1776 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 1777 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 1778 ierr = VecDestroy(&ts->Frhs);CHKERRQ(ierr); 1779 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1780 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 1781 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 1782 ierr = VecDestroyVecs(ts->nwork,&ts->work);CHKERRQ(ierr); 1783 1784 ts->setupcalled = PETSC_FALSE; 1785 PetscFunctionReturn(0); 1786 } 1787 1788 #undef __FUNCT__ 1789 #define __FUNCT__ "TSDestroy" 1790 /*@ 1791 TSDestroy - Destroys the timestepper context that was created 1792 with TSCreate(). 1793 1794 Collective on TS 1795 1796 Input Parameter: 1797 . ts - the TS context obtained from TSCreate() 1798 1799 Level: beginner 1800 1801 .keywords: TS, timestepper, destroy 1802 1803 .seealso: TSCreate(), TSSetUp(), TSSolve() 1804 @*/ 1805 PetscErrorCode TSDestroy(TS *ts) 1806 { 1807 PetscErrorCode ierr; 1808 1809 PetscFunctionBegin; 1810 if (!*ts) PetscFunctionReturn(0); 1811 PetscValidHeaderSpecific((*ts),TS_CLASSID,1); 1812 if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; PetscFunctionReturn(0);} 1813 1814 ierr = TSReset((*ts));CHKERRQ(ierr); 1815 1816 /* if memory was published with AMS then destroy it */ 1817 ierr = PetscObjectAMSViewOff((PetscObject)*ts);CHKERRQ(ierr); 1818 if ((*ts)->ops->destroy) {ierr = (*(*ts)->ops->destroy)((*ts));CHKERRQ(ierr);} 1819 1820 ierr = TSAdaptDestroy(&(*ts)->adapt);CHKERRQ(ierr); 1821 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 1822 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 1823 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 1824 1825 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 1826 PetscFunctionReturn(0); 1827 } 1828 1829 #undef __FUNCT__ 1830 #define __FUNCT__ "TSGetSNES" 1831 /*@ 1832 TSGetSNES - Returns the SNES (nonlinear solver) associated with 1833 a TS (timestepper) context. Valid only for nonlinear problems. 1834 1835 Not Collective, but SNES is parallel if TS is parallel 1836 1837 Input Parameter: 1838 . ts - the TS context obtained from TSCreate() 1839 1840 Output Parameter: 1841 . snes - the nonlinear solver context 1842 1843 Notes: 1844 The user can then directly manipulate the SNES context to set various 1845 options, etc. Likewise, the user can then extract and manipulate the 1846 KSP, KSP, and PC contexts as well. 1847 1848 TSGetSNES() does not work for integrators that do not use SNES; in 1849 this case TSGetSNES() returns NULL in snes. 1850 1851 Level: beginner 1852 1853 .keywords: timestep, get, SNES 1854 @*/ 1855 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 1856 { 1857 PetscErrorCode ierr; 1858 1859 PetscFunctionBegin; 1860 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1861 PetscValidPointer(snes,2); 1862 if (!ts->snes) { 1863 ierr = SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);CHKERRQ(ierr); 1864 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 1865 ierr = PetscLogObjectParent(ts,ts->snes);CHKERRQ(ierr); 1866 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 1867 if (ts->dm) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 1868 if (ts->problem_type == TS_LINEAR) { 1869 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 1870 } 1871 } 1872 *snes = ts->snes; 1873 PetscFunctionReturn(0); 1874 } 1875 1876 #undef __FUNCT__ 1877 #define __FUNCT__ "TSSetSNES" 1878 /*@ 1879 TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context 1880 1881 Collective 1882 1883 Input Parameter: 1884 + ts - the TS context obtained from TSCreate() 1885 - snes - the nonlinear solver context 1886 1887 Notes: 1888 Most users should have the TS created by calling TSGetSNES() 1889 1890 Level: developer 1891 1892 .keywords: timestep, set, SNES 1893 @*/ 1894 PetscErrorCode TSSetSNES(TS ts,SNES snes) 1895 { 1896 PetscErrorCode ierr; 1897 PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 1898 1899 PetscFunctionBegin; 1900 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1901 PetscValidHeaderSpecific(snes,SNES_CLASSID,2); 1902 ierr = PetscObjectReference((PetscObject)snes);CHKERRQ(ierr); 1903 ierr = SNESDestroy(&ts->snes);CHKERRQ(ierr); 1904 1905 ts->snes = snes; 1906 1907 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 1908 ierr = SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);CHKERRQ(ierr); 1909 if (func == SNESTSFormJacobian) { 1910 ierr = SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1911 } 1912 PetscFunctionReturn(0); 1913 } 1914 1915 #undef __FUNCT__ 1916 #define __FUNCT__ "TSGetKSP" 1917 /*@ 1918 TSGetKSP - Returns the KSP (linear solver) associated with 1919 a TS (timestepper) context. 1920 1921 Not Collective, but KSP is parallel if TS is parallel 1922 1923 Input Parameter: 1924 . ts - the TS context obtained from TSCreate() 1925 1926 Output Parameter: 1927 . ksp - the nonlinear solver context 1928 1929 Notes: 1930 The user can then directly manipulate the KSP context to set various 1931 options, etc. Likewise, the user can then extract and manipulate the 1932 KSP and PC contexts as well. 1933 1934 TSGetKSP() does not work for integrators that do not use KSP; 1935 in this case TSGetKSP() returns NULL in ksp. 1936 1937 Level: beginner 1938 1939 .keywords: timestep, get, KSP 1940 @*/ 1941 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 1942 { 1943 PetscErrorCode ierr; 1944 SNES snes; 1945 1946 PetscFunctionBegin; 1947 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1948 PetscValidPointer(ksp,2); 1949 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 1950 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 1951 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1952 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 1953 PetscFunctionReturn(0); 1954 } 1955 1956 /* ----------- Routines to set solver parameters ---------- */ 1957 1958 #undef __FUNCT__ 1959 #define __FUNCT__ "TSGetDuration" 1960 /*@ 1961 TSGetDuration - Gets the maximum number of timesteps to use and 1962 maximum time for iteration. 1963 1964 Not Collective 1965 1966 Input Parameters: 1967 + ts - the TS context obtained from TSCreate() 1968 . maxsteps - maximum number of iterations to use, or NULL 1969 - maxtime - final time to iterate to, or NULL 1970 1971 Level: intermediate 1972 1973 .keywords: TS, timestep, get, maximum, iterations, time 1974 @*/ 1975 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 1976 { 1977 PetscFunctionBegin; 1978 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1979 if (maxsteps) { 1980 PetscValidIntPointer(maxsteps,2); 1981 *maxsteps = ts->max_steps; 1982 } 1983 if (maxtime) { 1984 PetscValidScalarPointer(maxtime,3); 1985 *maxtime = ts->max_time; 1986 } 1987 PetscFunctionReturn(0); 1988 } 1989 1990 #undef __FUNCT__ 1991 #define __FUNCT__ "TSSetDuration" 1992 /*@ 1993 TSSetDuration - Sets the maximum number of timesteps to use and 1994 maximum time for iteration. 1995 1996 Logically Collective on TS 1997 1998 Input Parameters: 1999 + ts - the TS context obtained from TSCreate() 2000 . maxsteps - maximum number of iterations to use 2001 - maxtime - final time to iterate to 2002 2003 Options Database Keys: 2004 . -ts_max_steps <maxsteps> - Sets maxsteps 2005 . -ts_final_time <maxtime> - Sets maxtime 2006 2007 Notes: 2008 The default maximum number of iterations is 5000. Default time is 5.0 2009 2010 Level: intermediate 2011 2012 .keywords: TS, timestep, set, maximum, iterations 2013 2014 .seealso: TSSetExactFinalTime() 2015 @*/ 2016 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 2017 { 2018 PetscFunctionBegin; 2019 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2020 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 2021 PetscValidLogicalCollectiveReal(ts,maxtime,2); 2022 if (maxsteps >= 0) ts->max_steps = maxsteps; 2023 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 2024 PetscFunctionReturn(0); 2025 } 2026 2027 #undef __FUNCT__ 2028 #define __FUNCT__ "TSSetSolution" 2029 /*@ 2030 TSSetSolution - Sets the initial solution vector 2031 for use by the TS routines. 2032 2033 Logically Collective on TS and Vec 2034 2035 Input Parameters: 2036 + ts - the TS context obtained from TSCreate() 2037 - u - the solution vector 2038 2039 Level: beginner 2040 2041 .keywords: TS, timestep, set, solution, initial conditions 2042 @*/ 2043 PetscErrorCode TSSetSolution(TS ts,Vec u) 2044 { 2045 PetscErrorCode ierr; 2046 DM dm; 2047 2048 PetscFunctionBegin; 2049 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2050 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2051 ierr = PetscObjectReference((PetscObject)u);CHKERRQ(ierr); 2052 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 2053 2054 ts->vec_sol = u; 2055 2056 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2057 ierr = DMShellSetGlobalVector(dm,u);CHKERRQ(ierr); 2058 PetscFunctionReturn(0); 2059 } 2060 2061 #undef __FUNCT__ 2062 #define __FUNCT__ "TSSetPreStep" 2063 /*@C 2064 TSSetPreStep - Sets the general-purpose function 2065 called once at the beginning of each time step. 2066 2067 Logically Collective on TS 2068 2069 Input Parameters: 2070 + ts - The TS context obtained from TSCreate() 2071 - func - The function 2072 2073 Calling sequence of func: 2074 . func (TS ts); 2075 2076 Level: intermediate 2077 2078 Note: 2079 If a step is rejected, TSStep() will call this routine again before each attempt. 2080 The last completed time step number can be queried using TSGetTimeStepNumber(), the 2081 size of the step being attempted can be obtained using TSGetTimeStep(). 2082 2083 .keywords: TS, timestep 2084 .seealso: TSSetPreStage(), TSSetPostStage(), TSSetPostStep(), TSStep() 2085 @*/ 2086 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 2087 { 2088 PetscFunctionBegin; 2089 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2090 ts->prestep = func; 2091 PetscFunctionReturn(0); 2092 } 2093 2094 #undef __FUNCT__ 2095 #define __FUNCT__ "TSPreStep" 2096 /*@ 2097 TSPreStep - Runs the user-defined pre-step function. 2098 2099 Collective on TS 2100 2101 Input Parameters: 2102 . ts - The TS context obtained from TSCreate() 2103 2104 Notes: 2105 TSPreStep() is typically used within time stepping implementations, 2106 so most users would not generally call this routine themselves. 2107 2108 Level: developer 2109 2110 .keywords: TS, timestep 2111 .seealso: TSSetPreStep(), TSPreStage(), TSPostStage(), TSPostStep() 2112 @*/ 2113 PetscErrorCode TSPreStep(TS ts) 2114 { 2115 PetscErrorCode ierr; 2116 2117 PetscFunctionBegin; 2118 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2119 if (ts->prestep) { 2120 PetscStackCallStandard((*ts->prestep),(ts)); 2121 } 2122 PetscFunctionReturn(0); 2123 } 2124 2125 #undef __FUNCT__ 2126 #define __FUNCT__ "TSSetPreStage" 2127 /*@C 2128 TSSetPreStage - Sets the general-purpose function 2129 called once at the beginning of each stage. 2130 2131 Logically Collective on TS 2132 2133 Input Parameters: 2134 + ts - The TS context obtained from TSCreate() 2135 - func - The function 2136 2137 Calling sequence of func: 2138 . PetscErrorCode func(TS ts, PetscReal stagetime); 2139 2140 Level: intermediate 2141 2142 Note: 2143 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2144 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2145 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2146 2147 .keywords: TS, timestep 2148 .seealso: TSSetPostStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2149 @*/ 2150 PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal)) 2151 { 2152 PetscFunctionBegin; 2153 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2154 ts->prestage = func; 2155 PetscFunctionReturn(0); 2156 } 2157 2158 #undef __FUNCT__ 2159 #define __FUNCT__ "TSSetPostStage" 2160 /*@C 2161 TSSetPostStage - Sets the general-purpose function 2162 called once at the end of each stage. 2163 2164 Logically Collective on TS 2165 2166 Input Parameters: 2167 + ts - The TS context obtained from TSCreate() 2168 - func - The function 2169 2170 Calling sequence of func: 2171 . PetscErrorCode func(TS ts, PetscReal stagetime, PetscInt stageindex, Vec* Y); 2172 2173 Level: intermediate 2174 2175 Note: 2176 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2177 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2178 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2179 2180 .keywords: TS, timestep 2181 .seealso: TSSetPreStage(), TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2182 @*/ 2183 PetscErrorCode TSSetPostStage(TS ts, PetscErrorCode (*func)(TS,PetscReal,PetscInt,Vec*)) 2184 { 2185 PetscFunctionBegin; 2186 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2187 ts->poststage = func; 2188 PetscFunctionReturn(0); 2189 } 2190 2191 #undef __FUNCT__ 2192 #define __FUNCT__ "TSPreStage" 2193 /*@ 2194 TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage() 2195 2196 Collective on TS 2197 2198 Input Parameters: 2199 . ts - The TS context obtained from TSCreate() 2200 stagetime - The absolute time of the current stage 2201 2202 Notes: 2203 TSPreStage() is typically used within time stepping implementations, 2204 most users would not generally call this routine themselves. 2205 2206 Level: developer 2207 2208 .keywords: TS, timestep 2209 .seealso: TSPostStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 2210 @*/ 2211 PetscErrorCode TSPreStage(TS ts, PetscReal stagetime) 2212 { 2213 PetscErrorCode ierr; 2214 2215 PetscFunctionBegin; 2216 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2217 if (ts->prestage) { 2218 PetscStackCallStandard((*ts->prestage),(ts,stagetime)); 2219 } 2220 PetscFunctionReturn(0); 2221 } 2222 2223 #undef __FUNCT__ 2224 #define __FUNCT__ "TSPostStage" 2225 /*@ 2226 TSPostStage - Runs the user-defined post-stage function set using TSSetPostStage() 2227 2228 Collective on TS 2229 2230 Input Parameters: 2231 . ts - The TS context obtained from TSCreate() 2232 stagetime - The absolute time of the current stage 2233 stageindex - Stage number 2234 Y - Array of vectors (of size = total number 2235 of stages) with the stage solutions 2236 2237 Notes: 2238 TSPostStage() 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: TSPreStage(), TSSetPreStep(), TSPreStep(), TSPostStep() 2245 @*/ 2246 PetscErrorCode TSPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y) 2247 { 2248 PetscErrorCode ierr; 2249 2250 PetscFunctionBegin; 2251 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2252 if (ts->prestage) { 2253 PetscStackCallStandard((*ts->poststage),(ts,stagetime,stageindex,Y)); 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(), TSSetPostStage(), 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(), TSSetPostStage(), 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