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