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_AMS) 1280 #include <petscviewerams.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_AMS) 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_AMS) 1336 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERAMS,&isams);CHKERRQ(ierr); 1337 #endif 1338 if (iascii) { 1339 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer,"TS Object");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_AMS) 1395 } else if (isams) { 1396 if (((PetscObject)ts)->amsmem == -1) { 1397 ierr = PetscObjectViewAMS((PetscObject)ts,viewer);CHKERRQ(ierr); 1398 PetscStackCallAMS(AMS_Memory_take_access,(((PetscObject)ts)->amsmem)); 1399 PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)ts)->amsmem,"time step",&ts->steps,1,AMS_INT,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF)); 1400 PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)ts)->amsmem,"time",&ts->ptime,1,AMS_DOUBLE,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF)); 1401 PetscStackCallAMS(AMS_Memory_grant_access,(((PetscObject)ts)->amsmem)); 1402 } 1403 if (ts->ops->view) { 1404 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1405 } 1406 #endif 1407 } 1408 1409 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1410 ierr = PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);CHKERRQ(ierr); 1411 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1412 PetscFunctionReturn(0); 1413 } 1414 1415 1416 #undef __FUNCT__ 1417 #define __FUNCT__ "TSSetApplicationContext" 1418 /*@ 1419 TSSetApplicationContext - Sets an optional user-defined context for 1420 the timesteppers. 1421 1422 Logically Collective on TS 1423 1424 Input Parameters: 1425 + ts - the TS context obtained from TSCreate() 1426 - usrP - optional user context 1427 1428 Level: intermediate 1429 1430 .keywords: TS, timestep, set, application, context 1431 1432 .seealso: TSGetApplicationContext() 1433 @*/ 1434 PetscErrorCode TSSetApplicationContext(TS ts,void *usrP) 1435 { 1436 PetscFunctionBegin; 1437 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1438 ts->user = usrP; 1439 PetscFunctionReturn(0); 1440 } 1441 1442 #undef __FUNCT__ 1443 #define __FUNCT__ "TSGetApplicationContext" 1444 /*@ 1445 TSGetApplicationContext - Gets the user-defined context for the 1446 timestepper. 1447 1448 Not Collective 1449 1450 Input Parameter: 1451 . ts - the TS context obtained from TSCreate() 1452 1453 Output Parameter: 1454 . usrP - user context 1455 1456 Level: intermediate 1457 1458 .keywords: TS, timestep, get, application, context 1459 1460 .seealso: TSSetApplicationContext() 1461 @*/ 1462 PetscErrorCode TSGetApplicationContext(TS ts,void *usrP) 1463 { 1464 PetscFunctionBegin; 1465 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1466 *(void**)usrP = ts->user; 1467 PetscFunctionReturn(0); 1468 } 1469 1470 #undef __FUNCT__ 1471 #define __FUNCT__ "TSGetTimeStepNumber" 1472 /*@ 1473 TSGetTimeStepNumber - Gets the number of time steps completed. 1474 1475 Not Collective 1476 1477 Input Parameter: 1478 . ts - the TS context obtained from TSCreate() 1479 1480 Output Parameter: 1481 . iter - number of steps completed so far 1482 1483 Level: intermediate 1484 1485 .keywords: TS, timestep, get, iteration, number 1486 .seealso: TSGetTime(), TSGetTimeStep(), TSSetPreStep(), TSSetPreStage(), TSSetPostStep() 1487 @*/ 1488 PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt *iter) 1489 { 1490 PetscFunctionBegin; 1491 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1492 PetscValidIntPointer(iter,2); 1493 *iter = ts->steps; 1494 PetscFunctionReturn(0); 1495 } 1496 1497 #undef __FUNCT__ 1498 #define __FUNCT__ "TSSetInitialTimeStep" 1499 /*@ 1500 TSSetInitialTimeStep - Sets the initial timestep to be used, 1501 as well as the initial time. 1502 1503 Logically Collective on TS 1504 1505 Input Parameters: 1506 + ts - the TS context obtained from TSCreate() 1507 . initial_time - the initial time 1508 - time_step - the size of the timestep 1509 1510 Level: intermediate 1511 1512 .seealso: TSSetTimeStep(), TSGetTimeStep() 1513 1514 .keywords: TS, set, initial, timestep 1515 @*/ 1516 PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step) 1517 { 1518 PetscErrorCode ierr; 1519 1520 PetscFunctionBegin; 1521 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1522 ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr); 1523 ierr = TSSetTime(ts,initial_time);CHKERRQ(ierr); 1524 PetscFunctionReturn(0); 1525 } 1526 1527 #undef __FUNCT__ 1528 #define __FUNCT__ "TSSetTimeStep" 1529 /*@ 1530 TSSetTimeStep - Allows one to reset the timestep at any time, 1531 useful for simple pseudo-timestepping codes. 1532 1533 Logically Collective on TS 1534 1535 Input Parameters: 1536 + ts - the TS context obtained from TSCreate() 1537 - time_step - the size of the timestep 1538 1539 Level: intermediate 1540 1541 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1542 1543 .keywords: TS, set, timestep 1544 @*/ 1545 PetscErrorCode TSSetTimeStep(TS ts,PetscReal time_step) 1546 { 1547 PetscFunctionBegin; 1548 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1549 PetscValidLogicalCollectiveReal(ts,time_step,2); 1550 ts->time_step = time_step; 1551 ts->time_step_orig = time_step; 1552 PetscFunctionReturn(0); 1553 } 1554 1555 #undef __FUNCT__ 1556 #define __FUNCT__ "TSSetExactFinalTime" 1557 /*@ 1558 TSSetExactFinalTime - Determines whether to adapt the final time step to 1559 match the exact final time, interpolate solution to the exact final time, 1560 or just return at the final time TS computed. 1561 1562 Logically Collective on TS 1563 1564 Input Parameter: 1565 + ts - the time-step context 1566 - eftopt - exact final time option 1567 1568 Level: beginner 1569 1570 .seealso: TSExactFinalTimeOption 1571 @*/ 1572 PetscErrorCode TSSetExactFinalTime(TS ts,TSExactFinalTimeOption eftopt) 1573 { 1574 PetscFunctionBegin; 1575 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1576 PetscValidLogicalCollectiveEnum(ts,eftopt,2); 1577 ts->exact_final_time = eftopt; 1578 PetscFunctionReturn(0); 1579 } 1580 1581 #undef __FUNCT__ 1582 #define __FUNCT__ "TSGetTimeStep" 1583 /*@ 1584 TSGetTimeStep - Gets the current timestep size. 1585 1586 Not Collective 1587 1588 Input Parameter: 1589 . ts - the TS context obtained from TSCreate() 1590 1591 Output Parameter: 1592 . dt - the current timestep size 1593 1594 Level: intermediate 1595 1596 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1597 1598 .keywords: TS, get, timestep 1599 @*/ 1600 PetscErrorCode TSGetTimeStep(TS ts,PetscReal *dt) 1601 { 1602 PetscFunctionBegin; 1603 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1604 PetscValidRealPointer(dt,2); 1605 *dt = ts->time_step; 1606 PetscFunctionReturn(0); 1607 } 1608 1609 #undef __FUNCT__ 1610 #define __FUNCT__ "TSGetSolution" 1611 /*@ 1612 TSGetSolution - Returns the solution at the present timestep. It 1613 is valid to call this routine inside the function that you are evaluating 1614 in order to move to the new timestep. This vector not changed until 1615 the solution at the next timestep has been calculated. 1616 1617 Not Collective, but Vec returned is parallel if TS is parallel 1618 1619 Input Parameter: 1620 . ts - the TS context obtained from TSCreate() 1621 1622 Output Parameter: 1623 . v - the vector containing the solution 1624 1625 Level: intermediate 1626 1627 .seealso: TSGetTimeStep() 1628 1629 .keywords: TS, timestep, get, solution 1630 @*/ 1631 PetscErrorCode TSGetSolution(TS ts,Vec *v) 1632 { 1633 PetscFunctionBegin; 1634 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1635 PetscValidPointer(v,2); 1636 *v = ts->vec_sol; 1637 PetscFunctionReturn(0); 1638 } 1639 1640 /* ----- Routines to initialize and destroy a timestepper ---- */ 1641 #undef __FUNCT__ 1642 #define __FUNCT__ "TSSetProblemType" 1643 /*@ 1644 TSSetProblemType - Sets the type of problem to be solved. 1645 1646 Not collective 1647 1648 Input Parameters: 1649 + ts - The TS 1650 - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1651 .vb 1652 U_t - A U = 0 (linear) 1653 U_t - A(t) U = 0 (linear) 1654 F(t,U,U_t) = 0 (nonlinear) 1655 .ve 1656 1657 Level: beginner 1658 1659 .keywords: TS, problem type 1660 .seealso: TSSetUp(), TSProblemType, TS 1661 @*/ 1662 PetscErrorCode TSSetProblemType(TS ts, TSProblemType type) 1663 { 1664 PetscErrorCode ierr; 1665 1666 PetscFunctionBegin; 1667 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1668 ts->problem_type = type; 1669 if (type == TS_LINEAR) { 1670 SNES snes; 1671 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1672 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1673 } 1674 PetscFunctionReturn(0); 1675 } 1676 1677 #undef __FUNCT__ 1678 #define __FUNCT__ "TSGetProblemType" 1679 /*@C 1680 TSGetProblemType - Gets the type of problem to be solved. 1681 1682 Not collective 1683 1684 Input Parameter: 1685 . ts - The TS 1686 1687 Output Parameter: 1688 . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1689 .vb 1690 M U_t = A U 1691 M(t) U_t = A(t) U 1692 F(t,U,U_t) 1693 .ve 1694 1695 Level: beginner 1696 1697 .keywords: TS, problem type 1698 .seealso: TSSetUp(), TSProblemType, TS 1699 @*/ 1700 PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type) 1701 { 1702 PetscFunctionBegin; 1703 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1704 PetscValidIntPointer(type,2); 1705 *type = ts->problem_type; 1706 PetscFunctionReturn(0); 1707 } 1708 1709 #undef __FUNCT__ 1710 #define __FUNCT__ "TSSetUp" 1711 /*@ 1712 TSSetUp - Sets up the internal data structures for the later use 1713 of a timestepper. 1714 1715 Collective on TS 1716 1717 Input Parameter: 1718 . ts - the TS context obtained from TSCreate() 1719 1720 Notes: 1721 For basic use of the TS solvers the user need not explicitly call 1722 TSSetUp(), since these actions will automatically occur during 1723 the call to TSStep(). However, if one wishes to control this 1724 phase separately, TSSetUp() should be called after TSCreate() 1725 and optional routines of the form TSSetXXX(), but before TSStep(). 1726 1727 Level: advanced 1728 1729 .keywords: TS, timestep, setup 1730 1731 .seealso: TSCreate(), TSStep(), TSDestroy() 1732 @*/ 1733 PetscErrorCode TSSetUp(TS ts) 1734 { 1735 PetscErrorCode ierr; 1736 DM dm; 1737 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 1738 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 1739 TSIJacobian ijac; 1740 TSRHSJacobian rhsjac; 1741 1742 PetscFunctionBegin; 1743 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1744 if (ts->setupcalled) PetscFunctionReturn(0); 1745 1746 if (!((PetscObject)ts)->type_name) { 1747 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 1748 } 1749 1750 if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first"); 1751 1752 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 1753 1754 if (ts->rhsjacobian.reuse) { 1755 Mat Amat,Pmat; 1756 SNES snes; 1757 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1758 ierr = SNESGetJacobian(snes,&Amat,&Pmat,NULL,NULL);CHKERRQ(ierr); 1759 /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would 1760 * have displaced the RHS matrix */ 1761 if (Amat == ts->Arhs) { 1762 ierr = MatDuplicate(ts->Arhs,MAT_DO_NOT_COPY_VALUES,&Amat);CHKERRQ(ierr); 1763 ierr = SNESSetJacobian(snes,Amat,NULL,NULL,NULL);CHKERRQ(ierr); 1764 ierr = MatDestroy(&Amat);CHKERRQ(ierr); 1765 } 1766 if (Pmat == ts->Brhs) { 1767 ierr = MatDuplicate(ts->Brhs,MAT_DO_NOT_COPY_VALUES,&Pmat);CHKERRQ(ierr); 1768 ierr = SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);CHKERRQ(ierr); 1769 ierr = MatDestroy(&Pmat);CHKERRQ(ierr); 1770 } 1771 } 1772 1773 if (ts->ops->setup) { 1774 ierr = (*ts->ops->setup)(ts);CHKERRQ(ierr); 1775 } 1776 1777 /* in the case where we've set a DMTSFunction or what have you, we need the default SNESFunction 1778 to be set right but can't do it elsewhere due to the overreliance on ctx=ts. 1779 */ 1780 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1781 ierr = DMSNESGetFunction(dm,&func,NULL);CHKERRQ(ierr); 1782 if (!func) { 1783 ierr =DMSNESSetFunction(dm,SNESTSFormFunction,ts);CHKERRQ(ierr); 1784 } 1785 /* if the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it. 1786 Otherwise, the SNES will use coloring internally to form the Jacobian. 1787 */ 1788 ierr = DMSNESGetJacobian(dm,&jac,NULL);CHKERRQ(ierr); 1789 ierr = DMTSGetIJacobian(dm,&ijac,NULL);CHKERRQ(ierr); 1790 ierr = DMTSGetRHSJacobian(dm,&rhsjac,NULL);CHKERRQ(ierr); 1791 if (!jac && (ijac || rhsjac)) { 1792 ierr = DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1793 } 1794 ts->setupcalled = PETSC_TRUE; 1795 PetscFunctionReturn(0); 1796 } 1797 1798 #undef __FUNCT__ 1799 #define __FUNCT__ "TSReset" 1800 /*@ 1801 TSReset - Resets a TS context and removes any allocated Vecs and Mats. 1802 1803 Collective on TS 1804 1805 Input Parameter: 1806 . ts - the TS context obtained from TSCreate() 1807 1808 Level: beginner 1809 1810 .keywords: TS, timestep, reset 1811 1812 .seealso: TSCreate(), TSSetup(), TSDestroy() 1813 @*/ 1814 PetscErrorCode TSReset(TS ts) 1815 { 1816 PetscErrorCode ierr; 1817 1818 PetscFunctionBegin; 1819 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1820 if (ts->ops->reset) { 1821 ierr = (*ts->ops->reset)(ts);CHKERRQ(ierr); 1822 } 1823 if (ts->snes) {ierr = SNESReset(ts->snes);CHKERRQ(ierr);} 1824 1825 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 1826 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 1827 ierr = VecDestroy(&ts->Frhs);CHKERRQ(ierr); 1828 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1829 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 1830 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 1831 ierr = VecDestroyVecs(ts->nwork,&ts->work);CHKERRQ(ierr); 1832 1833 ts->setupcalled = PETSC_FALSE; 1834 PetscFunctionReturn(0); 1835 } 1836 1837 #undef __FUNCT__ 1838 #define __FUNCT__ "TSDestroy" 1839 /*@ 1840 TSDestroy - Destroys the timestepper context that was created 1841 with TSCreate(). 1842 1843 Collective on TS 1844 1845 Input Parameter: 1846 . ts - the TS context obtained from TSCreate() 1847 1848 Level: beginner 1849 1850 .keywords: TS, timestepper, destroy 1851 1852 .seealso: TSCreate(), TSSetUp(), TSSolve() 1853 @*/ 1854 PetscErrorCode TSDestroy(TS *ts) 1855 { 1856 PetscErrorCode ierr; 1857 1858 PetscFunctionBegin; 1859 if (!*ts) PetscFunctionReturn(0); 1860 PetscValidHeaderSpecific((*ts),TS_CLASSID,1); 1861 if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; PetscFunctionReturn(0);} 1862 1863 ierr = TSReset((*ts));CHKERRQ(ierr); 1864 1865 /* if memory was published with AMS then destroy it */ 1866 ierr = PetscObjectAMSViewOff((PetscObject)*ts);CHKERRQ(ierr); 1867 if ((*ts)->ops->destroy) {ierr = (*(*ts)->ops->destroy)((*ts));CHKERRQ(ierr);} 1868 1869 ierr = TSAdaptDestroy(&(*ts)->adapt);CHKERRQ(ierr); 1870 if ((*ts)->event) { 1871 ierr = TSEventMonitorDestroy(&(*ts)->event);CHKERRQ(ierr); 1872 } 1873 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 1874 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 1875 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 1876 1877 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 1878 PetscFunctionReturn(0); 1879 } 1880 1881 #undef __FUNCT__ 1882 #define __FUNCT__ "TSGetSNES" 1883 /*@ 1884 TSGetSNES - Returns the SNES (nonlinear solver) associated with 1885 a TS (timestepper) context. Valid only for nonlinear problems. 1886 1887 Not Collective, but SNES is parallel if TS is parallel 1888 1889 Input Parameter: 1890 . ts - the TS context obtained from TSCreate() 1891 1892 Output Parameter: 1893 . snes - the nonlinear solver context 1894 1895 Notes: 1896 The user can then directly manipulate the SNES context to set various 1897 options, etc. Likewise, the user can then extract and manipulate the 1898 KSP, KSP, and PC contexts as well. 1899 1900 TSGetSNES() does not work for integrators that do not use SNES; in 1901 this case TSGetSNES() returns NULL in snes. 1902 1903 Level: beginner 1904 1905 .keywords: timestep, get, SNES 1906 @*/ 1907 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 1908 { 1909 PetscErrorCode ierr; 1910 1911 PetscFunctionBegin; 1912 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1913 PetscValidPointer(snes,2); 1914 if (!ts->snes) { 1915 ierr = SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);CHKERRQ(ierr); 1916 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 1917 ierr = PetscLogObjectParent(ts,ts->snes);CHKERRQ(ierr); 1918 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 1919 if (ts->dm) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 1920 if (ts->problem_type == TS_LINEAR) { 1921 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 1922 } 1923 } 1924 *snes = ts->snes; 1925 PetscFunctionReturn(0); 1926 } 1927 1928 #undef __FUNCT__ 1929 #define __FUNCT__ "TSSetSNES" 1930 /*@ 1931 TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context 1932 1933 Collective 1934 1935 Input Parameter: 1936 + ts - the TS context obtained from TSCreate() 1937 - snes - the nonlinear solver context 1938 1939 Notes: 1940 Most users should have the TS created by calling TSGetSNES() 1941 1942 Level: developer 1943 1944 .keywords: timestep, set, SNES 1945 @*/ 1946 PetscErrorCode TSSetSNES(TS ts,SNES snes) 1947 { 1948 PetscErrorCode ierr; 1949 PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 1950 1951 PetscFunctionBegin; 1952 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1953 PetscValidHeaderSpecific(snes,SNES_CLASSID,2); 1954 ierr = PetscObjectReference((PetscObject)snes);CHKERRQ(ierr); 1955 ierr = SNESDestroy(&ts->snes);CHKERRQ(ierr); 1956 1957 ts->snes = snes; 1958 1959 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 1960 ierr = SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);CHKERRQ(ierr); 1961 if (func == SNESTSFormJacobian) { 1962 ierr = SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1963 } 1964 PetscFunctionReturn(0); 1965 } 1966 1967 #undef __FUNCT__ 1968 #define __FUNCT__ "TSGetKSP" 1969 /*@ 1970 TSGetKSP - Returns the KSP (linear solver) associated with 1971 a TS (timestepper) context. 1972 1973 Not Collective, but KSP is parallel if TS is parallel 1974 1975 Input Parameter: 1976 . ts - the TS context obtained from TSCreate() 1977 1978 Output Parameter: 1979 . ksp - the nonlinear solver context 1980 1981 Notes: 1982 The user can then directly manipulate the KSP context to set various 1983 options, etc. Likewise, the user can then extract and manipulate the 1984 KSP and PC contexts as well. 1985 1986 TSGetKSP() does not work for integrators that do not use KSP; 1987 in this case TSGetKSP() returns NULL in ksp. 1988 1989 Level: beginner 1990 1991 .keywords: timestep, get, KSP 1992 @*/ 1993 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 1994 { 1995 PetscErrorCode ierr; 1996 SNES snes; 1997 1998 PetscFunctionBegin; 1999 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2000 PetscValidPointer(ksp,2); 2001 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 2002 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 2003 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2004 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 2005 PetscFunctionReturn(0); 2006 } 2007 2008 /* ----------- Routines to set solver parameters ---------- */ 2009 2010 #undef __FUNCT__ 2011 #define __FUNCT__ "TSGetDuration" 2012 /*@ 2013 TSGetDuration - Gets the maximum number of timesteps to use and 2014 maximum time for iteration. 2015 2016 Not Collective 2017 2018 Input Parameters: 2019 + ts - the TS context obtained from TSCreate() 2020 . maxsteps - maximum number of iterations to use, or NULL 2021 - maxtime - final time to iterate to, or NULL 2022 2023 Level: intermediate 2024 2025 .keywords: TS, timestep, get, maximum, iterations, time 2026 @*/ 2027 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 2028 { 2029 PetscFunctionBegin; 2030 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2031 if (maxsteps) { 2032 PetscValidIntPointer(maxsteps,2); 2033 *maxsteps = ts->max_steps; 2034 } 2035 if (maxtime) { 2036 PetscValidScalarPointer(maxtime,3); 2037 *maxtime = ts->max_time; 2038 } 2039 PetscFunctionReturn(0); 2040 } 2041 2042 #undef __FUNCT__ 2043 #define __FUNCT__ "TSSetDuration" 2044 /*@ 2045 TSSetDuration - Sets the maximum number of timesteps to use and 2046 maximum time for iteration. 2047 2048 Logically Collective on TS 2049 2050 Input Parameters: 2051 + ts - the TS context obtained from TSCreate() 2052 . maxsteps - maximum number of iterations to use 2053 - maxtime - final time to iterate to 2054 2055 Options Database Keys: 2056 . -ts_max_steps <maxsteps> - Sets maxsteps 2057 . -ts_final_time <maxtime> - Sets maxtime 2058 2059 Notes: 2060 The default maximum number of iterations is 5000. Default time is 5.0 2061 2062 Level: intermediate 2063 2064 .keywords: TS, timestep, set, maximum, iterations 2065 2066 .seealso: TSSetExactFinalTime() 2067 @*/ 2068 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 2069 { 2070 PetscFunctionBegin; 2071 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2072 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 2073 PetscValidLogicalCollectiveReal(ts,maxtime,2); 2074 if (maxsteps >= 0) ts->max_steps = maxsteps; 2075 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 2076 PetscFunctionReturn(0); 2077 } 2078 2079 #undef __FUNCT__ 2080 #define __FUNCT__ "TSSetSolution" 2081 /*@ 2082 TSSetSolution - Sets the initial solution vector 2083 for use by the TS routines. 2084 2085 Logically Collective on TS and Vec 2086 2087 Input Parameters: 2088 + ts - the TS context obtained from TSCreate() 2089 - u - the solution vector 2090 2091 Level: beginner 2092 2093 .keywords: TS, timestep, set, solution, initial conditions 2094 @*/ 2095 PetscErrorCode TSSetSolution(TS ts,Vec u) 2096 { 2097 PetscErrorCode ierr; 2098 DM dm; 2099 2100 PetscFunctionBegin; 2101 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2102 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2103 ierr = PetscObjectReference((PetscObject)u);CHKERRQ(ierr); 2104 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 2105 2106 ts->vec_sol = u; 2107 2108 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2109 ierr = DMShellSetGlobalVector(dm,u);CHKERRQ(ierr); 2110 PetscFunctionReturn(0); 2111 } 2112 2113 #undef __FUNCT__ 2114 #define __FUNCT__ "TSSetPreStep" 2115 /*@C 2116 TSSetPreStep - Sets the general-purpose function 2117 called once at the beginning of each time step. 2118 2119 Logically Collective on TS 2120 2121 Input Parameters: 2122 + ts - The TS context obtained from TSCreate() 2123 - func - The function 2124 2125 Calling sequence of func: 2126 . func (TS ts); 2127 2128 Level: intermediate 2129 2130 Note: 2131 If a step is rejected, TSStep() will call this routine again before each attempt. 2132 The last completed time step number can be queried using TSGetTimeStepNumber(), the 2133 size of the step being attempted can be obtained using TSGetTimeStep(). 2134 2135 .keywords: TS, timestep 2136 .seealso: TSSetPreStage(), TSSetPostStep(), TSStep() 2137 @*/ 2138 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 2139 { 2140 PetscFunctionBegin; 2141 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2142 ts->prestep = func; 2143 PetscFunctionReturn(0); 2144 } 2145 2146 #undef __FUNCT__ 2147 #define __FUNCT__ "TSPreStep" 2148 /*@ 2149 TSPreStep - Runs the user-defined pre-step function. 2150 2151 Collective on TS 2152 2153 Input Parameters: 2154 . ts - The TS context obtained from TSCreate() 2155 2156 Notes: 2157 TSPreStep() is typically used within time stepping implementations, 2158 so most users would not generally call this routine themselves. 2159 2160 Level: developer 2161 2162 .keywords: TS, timestep 2163 .seealso: TSSetPreStep(), TSPreStage(), TSPostStep() 2164 @*/ 2165 PetscErrorCode TSPreStep(TS ts) 2166 { 2167 PetscErrorCode ierr; 2168 2169 PetscFunctionBegin; 2170 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2171 if (ts->prestep) { 2172 PetscStackCallStandard((*ts->prestep),(ts)); 2173 } 2174 PetscFunctionReturn(0); 2175 } 2176 2177 #undef __FUNCT__ 2178 #define __FUNCT__ "TSSetPreStage" 2179 /*@C 2180 TSSetPreStage - Sets the general-purpose function 2181 called once at the beginning of each stage. 2182 2183 Logically Collective on TS 2184 2185 Input Parameters: 2186 + ts - The TS context obtained from TSCreate() 2187 - func - The function 2188 2189 Calling sequence of func: 2190 . PetscErrorCode func(TS ts, PetscReal stagetime); 2191 2192 Level: intermediate 2193 2194 Note: 2195 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2196 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2197 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2198 2199 .keywords: TS, timestep 2200 .seealso: TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2201 @*/ 2202 PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal)) 2203 { 2204 PetscFunctionBegin; 2205 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2206 ts->prestage = func; 2207 PetscFunctionReturn(0); 2208 } 2209 2210 #undef __FUNCT__ 2211 #define __FUNCT__ "TSPreStage" 2212 /*@ 2213 TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage() 2214 2215 Collective on TS 2216 2217 Input Parameters: 2218 . ts - The TS context obtained from TSCreate() 2219 2220 Notes: 2221 TSPreStage() is typically used within time stepping implementations, 2222 most users would not generally call this routine themselves. 2223 2224 Level: developer 2225 2226 .keywords: TS, timestep 2227 .seealso: TSSetPreStep(), TSPreStep(), TSPostStep() 2228 @*/ 2229 PetscErrorCode TSPreStage(TS ts, PetscReal stagetime) 2230 { 2231 PetscErrorCode ierr; 2232 2233 PetscFunctionBegin; 2234 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2235 if (ts->prestage) { 2236 PetscStackCallStandard((*ts->prestage),(ts,stagetime)); 2237 } 2238 PetscFunctionReturn(0); 2239 } 2240 2241 #undef __FUNCT__ 2242 #define __FUNCT__ "TSSetPostStep" 2243 /*@C 2244 TSSetPostStep - Sets the general-purpose function 2245 called once at the end of each time step. 2246 2247 Logically Collective on TS 2248 2249 Input Parameters: 2250 + ts - The TS context obtained from TSCreate() 2251 - func - The function 2252 2253 Calling sequence of func: 2254 $ func (TS ts); 2255 2256 Level: intermediate 2257 2258 .keywords: TS, timestep 2259 .seealso: TSSetPreStep(), TSSetPreStage(), TSGetTimeStep(), TSGetTimeStepNumber(), TSGetTime() 2260 @*/ 2261 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 2262 { 2263 PetscFunctionBegin; 2264 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2265 ts->poststep = func; 2266 PetscFunctionReturn(0); 2267 } 2268 2269 #undef __FUNCT__ 2270 #define __FUNCT__ "TSPostStep" 2271 /*@ 2272 TSPostStep - Runs the user-defined post-step function. 2273 2274 Collective on TS 2275 2276 Input Parameters: 2277 . ts - The TS context obtained from TSCreate() 2278 2279 Notes: 2280 TSPostStep() is typically used within time stepping implementations, 2281 so most users would not generally call this routine themselves. 2282 2283 Level: developer 2284 2285 .keywords: TS, timestep 2286 @*/ 2287 PetscErrorCode TSPostStep(TS ts) 2288 { 2289 PetscErrorCode ierr; 2290 2291 PetscFunctionBegin; 2292 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2293 if (ts->poststep) { 2294 PetscStackCallStandard((*ts->poststep),(ts)); 2295 } 2296 PetscFunctionReturn(0); 2297 } 2298 2299 /* ------------ Routines to set performance monitoring options ----------- */ 2300 2301 #undef __FUNCT__ 2302 #define __FUNCT__ "TSMonitorSet" 2303 /*@C 2304 TSMonitorSet - Sets an ADDITIONAL function that is to be used at every 2305 timestep to display the iteration's progress. 2306 2307 Logically Collective on TS 2308 2309 Input Parameters: 2310 + ts - the TS context obtained from TSCreate() 2311 . monitor - monitoring routine 2312 . mctx - [optional] user-defined context for private data for the 2313 monitor routine (use NULL if no context is desired) 2314 - monitordestroy - [optional] routine that frees monitor context 2315 (may be NULL) 2316 2317 Calling sequence of monitor: 2318 $ int monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx) 2319 2320 + ts - the TS context 2321 . 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 2322 been interpolated to) 2323 . time - current time 2324 . u - current iterate 2325 - mctx - [optional] monitoring context 2326 2327 Notes: 2328 This routine adds an additional monitor to the list of monitors that 2329 already has been loaded. 2330 2331 Fortran notes: Only a single monitor function can be set for each TS object 2332 2333 Level: intermediate 2334 2335 .keywords: TS, timestep, set, monitor 2336 2337 .seealso: TSMonitorDefault(), TSMonitorCancel() 2338 @*/ 2339 PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**)) 2340 { 2341 PetscFunctionBegin; 2342 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2343 if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 2344 ts->monitor[ts->numbermonitors] = monitor; 2345 ts->monitordestroy[ts->numbermonitors] = mdestroy; 2346 ts->monitorcontext[ts->numbermonitors++] = (void*)mctx; 2347 PetscFunctionReturn(0); 2348 } 2349 2350 #undef __FUNCT__ 2351 #define __FUNCT__ "TSMonitorCancel" 2352 /*@C 2353 TSMonitorCancel - Clears all the monitors that have been set on a time-step object. 2354 2355 Logically Collective on TS 2356 2357 Input Parameters: 2358 . ts - the TS context obtained from TSCreate() 2359 2360 Notes: 2361 There is no way to remove a single, specific monitor. 2362 2363 Level: intermediate 2364 2365 .keywords: TS, timestep, set, monitor 2366 2367 .seealso: TSMonitorDefault(), TSMonitorSet() 2368 @*/ 2369 PetscErrorCode TSMonitorCancel(TS ts) 2370 { 2371 PetscErrorCode ierr; 2372 PetscInt i; 2373 2374 PetscFunctionBegin; 2375 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2376 for (i=0; i<ts->numbermonitors; i++) { 2377 if (ts->monitordestroy[i]) { 2378 ierr = (*ts->monitordestroy[i])(&ts->monitorcontext[i]);CHKERRQ(ierr); 2379 } 2380 } 2381 ts->numbermonitors = 0; 2382 PetscFunctionReturn(0); 2383 } 2384 2385 #undef __FUNCT__ 2386 #define __FUNCT__ "TSMonitorDefault" 2387 /*@ 2388 TSMonitorDefault - Sets the Default monitor 2389 2390 Level: intermediate 2391 2392 .keywords: TS, set, monitor 2393 2394 .seealso: TSMonitorDefault(), TSMonitorSet() 2395 @*/ 2396 PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy) 2397 { 2398 PetscErrorCode ierr; 2399 PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ts)); 2400 2401 PetscFunctionBegin; 2402 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2403 ierr = PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g\n",step,(double)ts->time_step,(double)ptime);CHKERRQ(ierr); 2404 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2405 PetscFunctionReturn(0); 2406 } 2407 2408 #undef __FUNCT__ 2409 #define __FUNCT__ "TSSetRetainStages" 2410 /*@ 2411 TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available. 2412 2413 Logically Collective on TS 2414 2415 Input Argument: 2416 . ts - time stepping context 2417 2418 Output Argument: 2419 . flg - PETSC_TRUE or PETSC_FALSE 2420 2421 Level: intermediate 2422 2423 .keywords: TS, set 2424 2425 .seealso: TSInterpolate(), TSSetPostStep() 2426 @*/ 2427 PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg) 2428 { 2429 PetscFunctionBegin; 2430 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2431 ts->retain_stages = flg; 2432 PetscFunctionReturn(0); 2433 } 2434 2435 #undef __FUNCT__ 2436 #define __FUNCT__ "TSInterpolate" 2437 /*@ 2438 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 2439 2440 Collective on TS 2441 2442 Input Argument: 2443 + ts - time stepping context 2444 - t - time to interpolate to 2445 2446 Output Argument: 2447 . U - state at given time 2448 2449 Notes: 2450 The user should call TSSetRetainStages() before taking a step in which interpolation will be requested. 2451 2452 Level: intermediate 2453 2454 Developer Notes: 2455 TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 2456 2457 .keywords: TS, set 2458 2459 .seealso: TSSetRetainStages(), TSSetPostStep() 2460 @*/ 2461 PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U) 2462 { 2463 PetscErrorCode ierr; 2464 2465 PetscFunctionBegin; 2466 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2467 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 2468 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); 2469 if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name); 2470 ierr = (*ts->ops->interpolate)(ts,t,U);CHKERRQ(ierr); 2471 PetscFunctionReturn(0); 2472 } 2473 2474 #undef __FUNCT__ 2475 #define __FUNCT__ "TSStep" 2476 /*@ 2477 TSStep - Steps one time step 2478 2479 Collective on TS 2480 2481 Input Parameter: 2482 . ts - the TS context obtained from TSCreate() 2483 2484 Level: intermediate 2485 2486 Notes: 2487 The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may 2488 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 2489 2490 This may over-step the final time provided in TSSetDuration() depending on the time-step used. TSSolve() interpolates to exactly the 2491 time provided in TSSetDuration(). One can use TSInterpolate() to determine an interpolated solution within the final timestep. 2492 2493 .keywords: TS, timestep, solve 2494 2495 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate() 2496 @*/ 2497 PetscErrorCode TSStep(TS ts) 2498 { 2499 PetscErrorCode ierr; 2500 2501 PetscFunctionBegin; 2502 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2503 ierr = TSSetUp(ts);CHKERRQ(ierr); 2504 2505 ts->reason = TS_CONVERGED_ITERATING; 2506 ts->ptime_prev = ts->ptime; 2507 2508 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2509 ierr = (*ts->ops->step)(ts);CHKERRQ(ierr); 2510 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2511 2512 ts->time_step_prev = ts->ptime - ts->ptime_prev; 2513 2514 if (ts->reason < 0) { 2515 if (ts->errorifstepfailed) { 2516 if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) { 2517 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]); 2518 } else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 2519 } 2520 } else if (!ts->reason) { 2521 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 2522 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 2523 } 2524 PetscFunctionReturn(0); 2525 } 2526 2527 #undef __FUNCT__ 2528 #define __FUNCT__ "TSEvaluateStep" 2529 /*@ 2530 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 2531 2532 Collective on TS 2533 2534 Input Arguments: 2535 + ts - time stepping context 2536 . order - desired order of accuracy 2537 - done - whether the step was evaluated at this order (pass NULL to generate an error if not available) 2538 2539 Output Arguments: 2540 . U - state at the end of the current step 2541 2542 Level: advanced 2543 2544 Notes: 2545 This function cannot be called until all stages have been evaluated. 2546 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. 2547 2548 .seealso: TSStep(), TSAdapt 2549 @*/ 2550 PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done) 2551 { 2552 PetscErrorCode ierr; 2553 2554 PetscFunctionBegin; 2555 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2556 PetscValidType(ts,1); 2557 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 2558 if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name); 2559 ierr = (*ts->ops->evaluatestep)(ts,order,U,done);CHKERRQ(ierr); 2560 PetscFunctionReturn(0); 2561 } 2562 2563 #undef __FUNCT__ 2564 #define __FUNCT__ "TSSolve" 2565 /*@ 2566 TSSolve - Steps the requested number of timesteps. 2567 2568 Collective on TS 2569 2570 Input Parameter: 2571 + ts - the TS context obtained from TSCreate() 2572 - u - the solution vector (can be null if TSSetSolution() was used, otherwise must contain the initial conditions) 2573 2574 Level: beginner 2575 2576 Notes: 2577 The final time returned by this function may be different from the time of the internally 2578 held state accessible by TSGetSolution() and TSGetTime() because the method may have 2579 stepped over the final time. 2580 2581 .keywords: TS, timestep, solve 2582 2583 .seealso: TSCreate(), TSSetSolution(), TSStep() 2584 @*/ 2585 PetscErrorCode TSSolve(TS ts,Vec u) 2586 { 2587 PetscBool flg; 2588 PetscViewer viewer; 2589 Vec solution; 2590 PetscErrorCode ierr; 2591 PetscViewerFormat format; 2592 2593 PetscFunctionBegin; 2594 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2595 if (u) PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2596 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 */ 2597 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2598 if (!ts->vec_sol || u == ts->vec_sol) { 2599 ierr = VecDuplicate(u,&solution);CHKERRQ(ierr); 2600 ierr = TSSetSolution(ts,solution);CHKERRQ(ierr); 2601 ierr = VecDestroy(&solution);CHKERRQ(ierr); /* grant ownership */ 2602 } 2603 ierr = VecCopy(u,ts->vec_sol);CHKERRQ(ierr); 2604 } else if (u) { 2605 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 2606 } 2607 ierr = TSSetUp(ts);CHKERRQ(ierr); 2608 /* reset time step and iteration counters */ 2609 ts->steps = 0; 2610 ts->ksp_its = 0; 2611 ts->snes_its = 0; 2612 ts->num_snes_failures = 0; 2613 ts->reject = 0; 2614 ts->reason = TS_CONVERGED_ITERATING; 2615 2616 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject)ts)->prefix,"-ts_view_pre",&viewer,&format,&flg);CHKERRQ(ierr); 2617 if (flg && !PetscPreLoadingOn) { 2618 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 2619 ierr = TSView(ts,viewer);CHKERRQ(ierr); 2620 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 2621 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2622 } 2623 2624 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 2625 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 2626 ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr); 2627 ts->solvetime = ts->ptime; 2628 } else { 2629 /* steps the requested number of timesteps. */ 2630 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 2631 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 2632 while (!ts->reason) { 2633 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2634 ierr = TSStep(ts);CHKERRQ(ierr); 2635 ierr = TSPostStep(ts);CHKERRQ(ierr); 2636 if (ts->event) { 2637 ierr = TSEventMonitor(ts);CHKERRQ(ierr); 2638 } 2639 } 2640 if (ts->exact_final_time == TS_EXACTFINALTIME_INTERPOLATE && ts->ptime > ts->max_time) { 2641 ierr = TSInterpolate(ts,ts->max_time,u);CHKERRQ(ierr); 2642 ts->solvetime = ts->max_time; 2643 solution = u; 2644 } else { 2645 if (u) {ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr);} 2646 ts->solvetime = ts->ptime; 2647 solution = ts->vec_sol; 2648 } 2649 ierr = TSMonitor(ts,ts->steps,ts->solvetime,solution);CHKERRQ(ierr); 2650 } 2651 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject)ts)->prefix,"-ts_view",&viewer,&format,&flg);CHKERRQ(ierr); 2652 if (flg && !PetscPreLoadingOn) { 2653 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 2654 ierr = TSView(ts,viewer);CHKERRQ(ierr); 2655 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 2656 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2657 } 2658 PetscFunctionReturn(0); 2659 } 2660 2661 #undef __FUNCT__ 2662 #define __FUNCT__ "TSMonitor" 2663 /*@ 2664 TSMonitor - Runs all user-provided monitor routines set using TSMonitorSet() 2665 2666 Collective on TS 2667 2668 Input Parameters: 2669 + ts - time stepping context obtained from TSCreate() 2670 . step - step number that has just completed 2671 . ptime - model time of the state 2672 - u - state at the current model time 2673 2674 Notes: 2675 TSMonitor() is typically used within the time stepping implementations. 2676 Users might call this function when using the TSStep() interface instead of TSSolve(). 2677 2678 Level: advanced 2679 2680 .keywords: TS, timestep 2681 @*/ 2682 PetscErrorCode TSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec u) 2683 { 2684 PetscErrorCode ierr; 2685 PetscInt i,n = ts->numbermonitors; 2686 2687 PetscFunctionBegin; 2688 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2689 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 2690 for (i=0; i<n; i++) { 2691 ierr = (*ts->monitor[i])(ts,step,ptime,u,ts->monitorcontext[i]);CHKERRQ(ierr); 2692 } 2693 PetscFunctionReturn(0); 2694 } 2695 2696 /* ------------------------------------------------------------------------*/ 2697 struct _n_TSMonitorLGCtx { 2698 PetscDrawLG lg; 2699 PetscInt howoften; /* when > 0 uses step % howoften, when negative only final solution plotted */ 2700 PetscInt ksp_its,snes_its; 2701 }; 2702 2703 2704 #undef __FUNCT__ 2705 #define __FUNCT__ "TSMonitorLGCtxCreate" 2706 /*@C 2707 TSMonitorLGCtxCreate - Creates a line graph context for use with 2708 TS to monitor the solution process graphically in various ways 2709 2710 Collective on TS 2711 2712 Input Parameters: 2713 + host - the X display to open, or null for the local machine 2714 . label - the title to put in the title bar 2715 . x, y - the screen coordinates of the upper left coordinate of the window 2716 . m, n - the screen width and height in pixels 2717 - howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time 2718 2719 Output Parameter: 2720 . ctx - the context 2721 2722 Options Database Key: 2723 + -ts_monitor_lg_timestep - automatically sets line graph monitor 2724 . -ts_monitor_lg_solution - 2725 . -ts_monitor_lg_error - 2726 . -ts_monitor_lg_ksp_iterations - 2727 . -ts_monitor_lg_snes_iterations - 2728 - -lg_indicate_data_points <true,false> - indicate the data points (at each time step) on the plot; default is true 2729 2730 Notes: 2731 Use TSMonitorLGCtxDestroy() to destroy. 2732 2733 Level: intermediate 2734 2735 .keywords: TS, monitor, line graph, residual, seealso 2736 2737 .seealso: TSMonitorLGTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError() 2738 2739 @*/ 2740 PetscErrorCode TSMonitorLGCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorLGCtx *ctx) 2741 { 2742 PetscDraw win; 2743 PetscErrorCode ierr; 2744 PetscBool flg = PETSC_TRUE; 2745 2746 PetscFunctionBegin; 2747 ierr = PetscNew(struct _n_TSMonitorLGCtx,ctx);CHKERRQ(ierr); 2748 ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr); 2749 ierr = PetscDrawSetFromOptions(win);CHKERRQ(ierr); 2750 ierr = PetscDrawLGCreate(win,1,&(*ctx)->lg);CHKERRQ(ierr); 2751 ierr = PetscOptionsGetBool(NULL,"-lg_indicate_data_points",&flg,NULL);CHKERRQ(ierr); 2752 if (flg) { 2753 ierr = PetscDrawLGIndicateDataPoints((*ctx)->lg);CHKERRQ(ierr); 2754 } 2755 ierr = PetscLogObjectParent((*ctx)->lg,win);CHKERRQ(ierr); 2756 (*ctx)->howoften = howoften; 2757 PetscFunctionReturn(0); 2758 } 2759 2760 #undef __FUNCT__ 2761 #define __FUNCT__ "TSMonitorLGTimeStep" 2762 PetscErrorCode TSMonitorLGTimeStep(TS ts,PetscInt step,PetscReal ptime,Vec v,void *monctx) 2763 { 2764 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 2765 PetscReal x = ptime,y; 2766 PetscErrorCode ierr; 2767 2768 PetscFunctionBegin; 2769 if (!step) { 2770 PetscDrawAxis axis; 2771 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 2772 ierr = PetscDrawAxisSetLabels(axis,"Timestep as function of time","Time","Time step");CHKERRQ(ierr); 2773 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 2774 } 2775 ierr = TSGetTimeStep(ts,&y);CHKERRQ(ierr); 2776 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 2777 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 2778 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 2779 } 2780 PetscFunctionReturn(0); 2781 } 2782 2783 #undef __FUNCT__ 2784 #define __FUNCT__ "TSMonitorLGCtxDestroy" 2785 /*@C 2786 TSMonitorLGCtxDestroy - Destroys a line graph context that was created 2787 with TSMonitorLGCtxCreate(). 2788 2789 Collective on TSMonitorLGCtx 2790 2791 Input Parameter: 2792 . ctx - the monitor context 2793 2794 Level: intermediate 2795 2796 .keywords: TS, monitor, line graph, destroy 2797 2798 .seealso: TSMonitorLGCtxCreate(), TSMonitorSet(), TSMonitorLGTimeStep(); 2799 @*/ 2800 PetscErrorCode TSMonitorLGCtxDestroy(TSMonitorLGCtx *ctx) 2801 { 2802 PetscDraw draw; 2803 PetscErrorCode ierr; 2804 2805 PetscFunctionBegin; 2806 ierr = PetscDrawLGGetDraw((*ctx)->lg,&draw);CHKERRQ(ierr); 2807 ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr); 2808 ierr = PetscDrawLGDestroy(&(*ctx)->lg);CHKERRQ(ierr); 2809 ierr = PetscFree(*ctx);CHKERRQ(ierr); 2810 PetscFunctionReturn(0); 2811 } 2812 2813 #undef __FUNCT__ 2814 #define __FUNCT__ "TSGetTime" 2815 /*@ 2816 TSGetTime - Gets the time of the most recently completed step. 2817 2818 Not Collective 2819 2820 Input Parameter: 2821 . ts - the TS context obtained from TSCreate() 2822 2823 Output Parameter: 2824 . t - the current time 2825 2826 Level: beginner 2827 2828 Note: 2829 When called during time step evaluation (e.g. during residual evaluation or via hooks set using TSSetPreStep(), 2830 TSSetPreStage(), or TSSetPostStep()), the time is the time at the start of the step being evaluated. 2831 2832 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 2833 2834 .keywords: TS, get, time 2835 @*/ 2836 PetscErrorCode TSGetTime(TS ts,PetscReal *t) 2837 { 2838 PetscFunctionBegin; 2839 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2840 PetscValidRealPointer(t,2); 2841 *t = ts->ptime; 2842 PetscFunctionReturn(0); 2843 } 2844 2845 #undef __FUNCT__ 2846 #define __FUNCT__ "TSSetTime" 2847 /*@ 2848 TSSetTime - Allows one to reset the time. 2849 2850 Logically Collective on TS 2851 2852 Input Parameters: 2853 + ts - the TS context obtained from TSCreate() 2854 - time - the time 2855 2856 Level: intermediate 2857 2858 .seealso: TSGetTime(), TSSetDuration() 2859 2860 .keywords: TS, set, time 2861 @*/ 2862 PetscErrorCode TSSetTime(TS ts, PetscReal t) 2863 { 2864 PetscFunctionBegin; 2865 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2866 PetscValidLogicalCollectiveReal(ts,t,2); 2867 ts->ptime = t; 2868 PetscFunctionReturn(0); 2869 } 2870 2871 #undef __FUNCT__ 2872 #define __FUNCT__ "TSSetOptionsPrefix" 2873 /*@C 2874 TSSetOptionsPrefix - Sets the prefix used for searching for all 2875 TS options in the database. 2876 2877 Logically Collective on TS 2878 2879 Input Parameter: 2880 + ts - The TS context 2881 - prefix - The prefix to prepend to all option names 2882 2883 Notes: 2884 A hyphen (-) must NOT be given at the beginning of the prefix name. 2885 The first character of all runtime options is AUTOMATICALLY the 2886 hyphen. 2887 2888 Level: advanced 2889 2890 .keywords: TS, set, options, prefix, database 2891 2892 .seealso: TSSetFromOptions() 2893 2894 @*/ 2895 PetscErrorCode TSSetOptionsPrefix(TS ts,const char prefix[]) 2896 { 2897 PetscErrorCode ierr; 2898 SNES snes; 2899 2900 PetscFunctionBegin; 2901 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2902 ierr = PetscObjectSetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2903 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2904 ierr = SNESSetOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2905 PetscFunctionReturn(0); 2906 } 2907 2908 2909 #undef __FUNCT__ 2910 #define __FUNCT__ "TSAppendOptionsPrefix" 2911 /*@C 2912 TSAppendOptionsPrefix - Appends to the prefix used for searching for all 2913 TS options in the database. 2914 2915 Logically Collective on TS 2916 2917 Input Parameter: 2918 + ts - The TS context 2919 - prefix - The prefix to prepend to all option names 2920 2921 Notes: 2922 A hyphen (-) must NOT be given at the beginning of the prefix name. 2923 The first character of all runtime options is AUTOMATICALLY the 2924 hyphen. 2925 2926 Level: advanced 2927 2928 .keywords: TS, append, options, prefix, database 2929 2930 .seealso: TSGetOptionsPrefix() 2931 2932 @*/ 2933 PetscErrorCode TSAppendOptionsPrefix(TS ts,const char prefix[]) 2934 { 2935 PetscErrorCode ierr; 2936 SNES snes; 2937 2938 PetscFunctionBegin; 2939 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2940 ierr = PetscObjectAppendOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2941 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2942 ierr = SNESAppendOptionsPrefix(snes,prefix);CHKERRQ(ierr); 2943 PetscFunctionReturn(0); 2944 } 2945 2946 #undef __FUNCT__ 2947 #define __FUNCT__ "TSGetOptionsPrefix" 2948 /*@C 2949 TSGetOptionsPrefix - Sets the prefix used for searching for all 2950 TS options in the database. 2951 2952 Not Collective 2953 2954 Input Parameter: 2955 . ts - The TS context 2956 2957 Output Parameter: 2958 . prefix - A pointer to the prefix string used 2959 2960 Notes: On the fortran side, the user should pass in a string 'prifix' of 2961 sufficient length to hold the prefix. 2962 2963 Level: intermediate 2964 2965 .keywords: TS, get, options, prefix, database 2966 2967 .seealso: TSAppendOptionsPrefix() 2968 @*/ 2969 PetscErrorCode TSGetOptionsPrefix(TS ts,const char *prefix[]) 2970 { 2971 PetscErrorCode ierr; 2972 2973 PetscFunctionBegin; 2974 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2975 PetscValidPointer(prefix,2); 2976 ierr = PetscObjectGetOptionsPrefix((PetscObject)ts,prefix);CHKERRQ(ierr); 2977 PetscFunctionReturn(0); 2978 } 2979 2980 #undef __FUNCT__ 2981 #define __FUNCT__ "TSGetRHSJacobian" 2982 /*@C 2983 TSGetRHSJacobian - Returns the Jacobian J at the present timestep. 2984 2985 Not Collective, but parallel objects are returned if TS is parallel 2986 2987 Input Parameter: 2988 . ts - The TS context obtained from TSCreate() 2989 2990 Output Parameters: 2991 + Amat - The (approximate) Jacobian J of G, where U_t = G(U,t) (or NULL) 2992 . Pmat - The matrix from which the preconditioner is constructed, usually the same as Amat (or NULL) 2993 . func - Function to compute the Jacobian of the RHS (or NULL) 2994 - ctx - User-defined context for Jacobian evaluation routine (or NULL) 2995 2996 Notes: You can pass in NULL for any return argument you do not need. 2997 2998 Level: intermediate 2999 3000 .seealso: TSGetTimeStep(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3001 3002 .keywords: TS, timestep, get, matrix, Jacobian 3003 @*/ 3004 PetscErrorCode TSGetRHSJacobian(TS ts,Mat *Amat,Mat *Pmat,TSRHSJacobian *func,void **ctx) 3005 { 3006 PetscErrorCode ierr; 3007 SNES snes; 3008 DM dm; 3009 3010 PetscFunctionBegin; 3011 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3012 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3013 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3014 ierr = DMTSGetRHSJacobian(dm,func,ctx);CHKERRQ(ierr); 3015 PetscFunctionReturn(0); 3016 } 3017 3018 #undef __FUNCT__ 3019 #define __FUNCT__ "TSGetIJacobian" 3020 /*@C 3021 TSGetIJacobian - Returns the implicit Jacobian at the present timestep. 3022 3023 Not Collective, but parallel objects are returned if TS is parallel 3024 3025 Input Parameter: 3026 . ts - The TS context obtained from TSCreate() 3027 3028 Output Parameters: 3029 + Amat - The (approximate) Jacobian of F(t,U,U_t) 3030 . Pmat - The matrix from which the preconditioner is constructed, often the same as Amat 3031 . f - The function to compute the matrices 3032 - ctx - User-defined context for Jacobian evaluation routine 3033 3034 Notes: You can pass in NULL for any return argument you do not need. 3035 3036 Level: advanced 3037 3038 .seealso: TSGetTimeStep(), TSGetRHSJacobian(), TSGetMatrices(), TSGetTime(), TSGetTimeStepNumber() 3039 3040 .keywords: TS, timestep, get, matrix, Jacobian 3041 @*/ 3042 PetscErrorCode TSGetIJacobian(TS ts,Mat *Amat,Mat *Pmat,TSIJacobian *f,void **ctx) 3043 { 3044 PetscErrorCode ierr; 3045 SNES snes; 3046 DM dm; 3047 3048 PetscFunctionBegin; 3049 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3050 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 3051 ierr = SNESGetJacobian(snes,Amat,Pmat,NULL,NULL);CHKERRQ(ierr); 3052 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 3053 ierr = DMTSGetIJacobian(dm,f,ctx);CHKERRQ(ierr); 3054 PetscFunctionReturn(0); 3055 } 3056 3057 3058 #undef __FUNCT__ 3059 #define __FUNCT__ "TSMonitorDrawSolution" 3060 /*@C 3061 TSMonitorDrawSolution - Monitors progress of the TS solvers by calling 3062 VecView() for the solution at each timestep 3063 3064 Collective on TS 3065 3066 Input Parameters: 3067 + ts - the TS context 3068 . step - current time-step 3069 . ptime - current time 3070 - dummy - either a viewer or NULL 3071 3072 Options Database: 3073 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3074 3075 Notes: the initial solution and current solution are not displayed with a common axis scaling so generally the option -ts_monitor_draw_solution_initial 3076 will look bad 3077 3078 Level: intermediate 3079 3080 .keywords: TS, vector, monitor, view 3081 3082 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3083 @*/ 3084 PetscErrorCode TSMonitorDrawSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3085 { 3086 PetscErrorCode ierr; 3087 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3088 PetscDraw draw; 3089 3090 PetscFunctionBegin; 3091 if (!step && ictx->showinitial) { 3092 if (!ictx->initialsolution) { 3093 ierr = VecDuplicate(u,&ictx->initialsolution);CHKERRQ(ierr); 3094 } 3095 ierr = VecCopy(u,ictx->initialsolution);CHKERRQ(ierr); 3096 } 3097 if (!(((ictx->howoften > 0) && (!(step % ictx->howoften))) || ((ictx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 3098 3099 if (ictx->showinitial) { 3100 PetscReal pause; 3101 ierr = PetscViewerDrawGetPause(ictx->viewer,&pause);CHKERRQ(ierr); 3102 ierr = PetscViewerDrawSetPause(ictx->viewer,0.0);CHKERRQ(ierr); 3103 ierr = VecView(ictx->initialsolution,ictx->viewer);CHKERRQ(ierr); 3104 ierr = PetscViewerDrawSetPause(ictx->viewer,pause);CHKERRQ(ierr); 3105 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_TRUE);CHKERRQ(ierr); 3106 } 3107 ierr = VecView(u,ictx->viewer);CHKERRQ(ierr); 3108 if (ictx->showtimestepandtime) { 3109 PetscReal xl,yl,xr,yr,tw,w,h; 3110 char time[32]; 3111 size_t len; 3112 3113 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3114 ierr = PetscSNPrintf(time,32,"Timestep %d Time %f",(int)step,(double)ptime);CHKERRQ(ierr); 3115 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3116 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3117 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3118 w = xl + .5*(xr - xl) - .5*len*tw; 3119 h = yl + .95*(yr - yl); 3120 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3121 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3122 } 3123 3124 if (ictx->showinitial) { 3125 ierr = PetscViewerDrawSetHold(ictx->viewer,PETSC_FALSE);CHKERRQ(ierr); 3126 } 3127 PetscFunctionReturn(0); 3128 } 3129 3130 #undef __FUNCT__ 3131 #define __FUNCT__ "TSMonitorDrawSolutionPhase" 3132 /*@C 3133 TSMonitorDrawSolutionPhase - Monitors progress of the TS solvers by plotting the solution as a phase diagram 3134 3135 Collective on TS 3136 3137 Input Parameters: 3138 + ts - the TS context 3139 . step - current time-step 3140 . ptime - current time 3141 - dummy - either a viewer or NULL 3142 3143 Level: intermediate 3144 3145 .keywords: TS, vector, monitor, view 3146 3147 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3148 @*/ 3149 PetscErrorCode TSMonitorDrawSolutionPhase(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3150 { 3151 PetscErrorCode ierr; 3152 TSMonitorDrawCtx ictx = (TSMonitorDrawCtx)dummy; 3153 PetscDraw draw; 3154 MPI_Comm comm; 3155 PetscInt n; 3156 PetscMPIInt size; 3157 PetscReal xl,yl,xr,yr,tw,w,h; 3158 char time[32]; 3159 size_t len; 3160 const PetscScalar *U; 3161 3162 PetscFunctionBegin; 3163 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 3164 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3165 if (size != 1) SETERRQ(comm,PETSC_ERR_SUP,"Only allowed for sequential runs"); 3166 ierr = VecGetSize(u,&n);CHKERRQ(ierr); 3167 if (n != 2) SETERRQ(comm,PETSC_ERR_SUP,"Only for ODEs with two unknowns"); 3168 3169 ierr = PetscViewerDrawGetDraw(ictx->viewer,0,&draw);CHKERRQ(ierr); 3170 3171 ierr = VecGetArrayRead(u,&U);CHKERRQ(ierr); 3172 ierr = PetscDrawAxisGetLimits(ictx->axis,&xl,&xr,&yl,&yr);CHKERRQ(ierr); 3173 if ((PetscRealPart(U[0]) < xl) || (PetscRealPart(U[1]) < yl) || (PetscRealPart(U[0]) > xr) || (PetscRealPart(U[1]) > yr)) { 3174 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3175 PetscFunctionReturn(0); 3176 } 3177 if (!step) ictx->color++; 3178 ierr = PetscDrawPoint(draw,PetscRealPart(U[0]),PetscRealPart(U[1]),ictx->color);CHKERRQ(ierr); 3179 ierr = VecRestoreArrayRead(u,&U);CHKERRQ(ierr); 3180 3181 if (ictx->showtimestepandtime) { 3182 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 3183 ierr = PetscSNPrintf(time,32,"Timestep %d Time %f",(int)step,(double)ptime);CHKERRQ(ierr); 3184 ierr = PetscStrlen(time,&len);CHKERRQ(ierr); 3185 ierr = PetscDrawStringGetSize(draw,&tw,NULL);CHKERRQ(ierr); 3186 w = xl + .5*(xr - xl) - .5*len*tw; 3187 h = yl + .95*(yr - yl); 3188 ierr = PetscDrawString(draw,w,h,PETSC_DRAW_BLACK,time);CHKERRQ(ierr); 3189 } 3190 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 3191 PetscFunctionReturn(0); 3192 } 3193 3194 3195 #undef __FUNCT__ 3196 #define __FUNCT__ "TSMonitorDrawCtxDestroy" 3197 /*@C 3198 TSMonitorDrawCtxDestroy - Destroys the monitor context for TSMonitorDrawSolution() 3199 3200 Collective on TS 3201 3202 Input Parameters: 3203 . ctx - the monitor context 3204 3205 Level: intermediate 3206 3207 .keywords: TS, vector, monitor, view 3208 3209 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawSolution(), TSMonitorDrawError() 3210 @*/ 3211 PetscErrorCode TSMonitorDrawCtxDestroy(TSMonitorDrawCtx *ictx) 3212 { 3213 PetscErrorCode ierr; 3214 3215 PetscFunctionBegin; 3216 ierr = PetscDrawAxisDestroy(&(*ictx)->axis);CHKERRQ(ierr); 3217 ierr = PetscViewerDestroy(&(*ictx)->viewer);CHKERRQ(ierr); 3218 ierr = VecDestroy(&(*ictx)->initialsolution);CHKERRQ(ierr); 3219 ierr = PetscFree(*ictx);CHKERRQ(ierr); 3220 PetscFunctionReturn(0); 3221 } 3222 3223 #undef __FUNCT__ 3224 #define __FUNCT__ "TSMonitorDrawCtxCreate" 3225 /*@C 3226 TSMonitorDrawCtxCreate - Creates the monitor context for TSMonitorDrawCtx 3227 3228 Collective on TS 3229 3230 Input Parameter: 3231 . ts - time-step context 3232 3233 Output Patameter: 3234 . ctx - the monitor context 3235 3236 Options Database: 3237 . -ts_monitor_draw_solution_initial - show initial solution as well as current solution 3238 3239 Level: intermediate 3240 3241 .keywords: TS, vector, monitor, view 3242 3243 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSMonitorDrawCtx() 3244 @*/ 3245 PetscErrorCode TSMonitorDrawCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorDrawCtx *ctx) 3246 { 3247 PetscErrorCode ierr; 3248 3249 PetscFunctionBegin; 3250 ierr = PetscNew(struct _n_TSMonitorDrawCtx,ctx);CHKERRQ(ierr); 3251 ierr = PetscViewerDrawOpen(comm,host,label,x,y,m,n,&(*ctx)->viewer);CHKERRQ(ierr); 3252 ierr = PetscViewerSetFromOptions((*ctx)->viewer);CHKERRQ(ierr); 3253 3254 (*ctx)->howoften = howoften; 3255 (*ctx)->showinitial = PETSC_FALSE; 3256 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_initial",&(*ctx)->showinitial,NULL);CHKERRQ(ierr); 3257 3258 (*ctx)->showtimestepandtime = PETSC_FALSE; 3259 ierr = PetscOptionsGetBool(NULL,"-ts_monitor_draw_solution_show_time",&(*ctx)->showtimestepandtime,NULL);CHKERRQ(ierr); 3260 (*ctx)->color = PETSC_DRAW_WHITE; 3261 PetscFunctionReturn(0); 3262 } 3263 3264 #undef __FUNCT__ 3265 #define __FUNCT__ "TSMonitorDrawError" 3266 /*@C 3267 TSMonitorDrawError - Monitors progress of the TS solvers by calling 3268 VecView() for the error at each timestep 3269 3270 Collective on TS 3271 3272 Input Parameters: 3273 + ts - the TS context 3274 . step - current time-step 3275 . ptime - current time 3276 - dummy - either a viewer or NULL 3277 3278 Level: intermediate 3279 3280 .keywords: TS, vector, monitor, view 3281 3282 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3283 @*/ 3284 PetscErrorCode TSMonitorDrawError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 3285 { 3286 PetscErrorCode ierr; 3287 TSMonitorDrawCtx ctx = (TSMonitorDrawCtx)dummy; 3288 PetscViewer viewer = ctx->viewer; 3289 Vec work; 3290 3291 PetscFunctionBegin; 3292 if (!(((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason))) PetscFunctionReturn(0); 3293 ierr = VecDuplicate(u,&work);CHKERRQ(ierr); 3294 ierr = TSComputeSolutionFunction(ts,ptime,work);CHKERRQ(ierr); 3295 ierr = VecAXPY(work,-1.0,u);CHKERRQ(ierr); 3296 ierr = VecView(work,viewer);CHKERRQ(ierr); 3297 ierr = VecDestroy(&work);CHKERRQ(ierr); 3298 PetscFunctionReturn(0); 3299 } 3300 3301 #include <petsc-private/dmimpl.h> 3302 #undef __FUNCT__ 3303 #define __FUNCT__ "TSSetDM" 3304 /*@ 3305 TSSetDM - Sets the DM that may be used by some preconditioners 3306 3307 Logically Collective on TS and DM 3308 3309 Input Parameters: 3310 + ts - the preconditioner context 3311 - dm - the dm 3312 3313 Level: intermediate 3314 3315 3316 .seealso: TSGetDM(), SNESSetDM(), SNESGetDM() 3317 @*/ 3318 PetscErrorCode TSSetDM(TS ts,DM dm) 3319 { 3320 PetscErrorCode ierr; 3321 SNES snes; 3322 DMTS tsdm; 3323 3324 PetscFunctionBegin; 3325 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3326 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 3327 if (ts->dm) { /* Move the DMTS context over to the new DM unless the new DM already has one */ 3328 if (ts->dm->dmts && !dm->dmts) { 3329 ierr = DMCopyDMTS(ts->dm,dm);CHKERRQ(ierr); 3330 ierr = DMGetDMTS(ts->dm,&tsdm);CHKERRQ(ierr); 3331 if (tsdm->originaldm == ts->dm) { /* Grant write privileges to the replacement DM */ 3332 tsdm->originaldm = dm; 3333 } 3334 } 3335 ierr = DMDestroy(&ts->dm);CHKERRQ(ierr); 3336 } 3337 ts->dm = dm; 3338 3339 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 3340 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 3341 PetscFunctionReturn(0); 3342 } 3343 3344 #undef __FUNCT__ 3345 #define __FUNCT__ "TSGetDM" 3346 /*@ 3347 TSGetDM - Gets the DM that may be used by some preconditioners 3348 3349 Not Collective 3350 3351 Input Parameter: 3352 . ts - the preconditioner context 3353 3354 Output Parameter: 3355 . dm - the dm 3356 3357 Level: intermediate 3358 3359 3360 .seealso: TSSetDM(), SNESSetDM(), SNESGetDM() 3361 @*/ 3362 PetscErrorCode TSGetDM(TS ts,DM *dm) 3363 { 3364 PetscErrorCode ierr; 3365 3366 PetscFunctionBegin; 3367 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3368 if (!ts->dm) { 3369 ierr = DMShellCreate(PetscObjectComm((PetscObject)ts),&ts->dm);CHKERRQ(ierr); 3370 if (ts->snes) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 3371 } 3372 *dm = ts->dm; 3373 PetscFunctionReturn(0); 3374 } 3375 3376 #undef __FUNCT__ 3377 #define __FUNCT__ "SNESTSFormFunction" 3378 /*@ 3379 SNESTSFormFunction - Function to evaluate nonlinear residual 3380 3381 Logically Collective on SNES 3382 3383 Input Parameter: 3384 + snes - nonlinear solver 3385 . U - the current state at which to evaluate the residual 3386 - ctx - user context, must be a TS 3387 3388 Output Parameter: 3389 . F - the nonlinear residual 3390 3391 Notes: 3392 This function is not normally called by users and is automatically registered with the SNES used by TS. 3393 It is most frequently passed to MatFDColoringSetFunction(). 3394 3395 Level: advanced 3396 3397 .seealso: SNESSetFunction(), MatFDColoringSetFunction() 3398 @*/ 3399 PetscErrorCode SNESTSFormFunction(SNES snes,Vec U,Vec F,void *ctx) 3400 { 3401 TS ts = (TS)ctx; 3402 PetscErrorCode ierr; 3403 3404 PetscFunctionBegin; 3405 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3406 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3407 PetscValidHeaderSpecific(F,VEC_CLASSID,3); 3408 PetscValidHeaderSpecific(ts,TS_CLASSID,4); 3409 ierr = (ts->ops->snesfunction)(snes,U,F,ts);CHKERRQ(ierr); 3410 PetscFunctionReturn(0); 3411 } 3412 3413 #undef __FUNCT__ 3414 #define __FUNCT__ "SNESTSFormJacobian" 3415 /*@ 3416 SNESTSFormJacobian - Function to evaluate the Jacobian 3417 3418 Collective on SNES 3419 3420 Input Parameter: 3421 + snes - nonlinear solver 3422 . U - the current state at which to evaluate the residual 3423 - ctx - user context, must be a TS 3424 3425 Output Parameter: 3426 + A - the Jacobian 3427 . B - the preconditioning matrix (may be the same as A) 3428 - flag - indicates any structure change in the matrix 3429 3430 Notes: 3431 This function is not normally called by users and is automatically registered with the SNES used by TS. 3432 3433 Level: developer 3434 3435 .seealso: SNESSetJacobian() 3436 @*/ 3437 PetscErrorCode SNESTSFormJacobian(SNES snes,Vec U,Mat *A,Mat *B,MatStructure *flag,void *ctx) 3438 { 3439 TS ts = (TS)ctx; 3440 PetscErrorCode ierr; 3441 3442 PetscFunctionBegin; 3443 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3444 PetscValidHeaderSpecific(U,VEC_CLASSID,2); 3445 PetscValidPointer(A,3); 3446 PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 3447 PetscValidPointer(B,4); 3448 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); 3449 PetscValidPointer(flag,5); 3450 PetscValidHeaderSpecific(ts,TS_CLASSID,6); 3451 ierr = (ts->ops->snesjacobian)(snes,U,A,B,flag,ts);CHKERRQ(ierr); 3452 PetscFunctionReturn(0); 3453 } 3454 3455 #undef __FUNCT__ 3456 #define __FUNCT__ "TSComputeRHSFunctionLinear" 3457 /*@C 3458 TSComputeRHSFunctionLinear - Evaluate the right hand side via the user-provided Jacobian, for linear problems only 3459 3460 Collective on TS 3461 3462 Input Arguments: 3463 + ts - time stepping context 3464 . t - time at which to evaluate 3465 . U - state at which to evaluate 3466 - ctx - context 3467 3468 Output Arguments: 3469 . F - right hand side 3470 3471 Level: intermediate 3472 3473 Notes: 3474 This function is intended to be passed to TSSetRHSFunction() to evaluate the right hand side for linear problems. 3475 The matrix (and optionally the evaluation context) should be passed to TSSetRHSJacobian(). 3476 3477 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 3478 @*/ 3479 PetscErrorCode TSComputeRHSFunctionLinear(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 3480 { 3481 PetscErrorCode ierr; 3482 Mat Arhs,Brhs; 3483 MatStructure flg2; 3484 3485 PetscFunctionBegin; 3486 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 3487 ierr = TSComputeRHSJacobian(ts,t,U,&Arhs,&Brhs,&flg2);CHKERRQ(ierr); 3488 ierr = MatMult(Arhs,U,F);CHKERRQ(ierr); 3489 PetscFunctionReturn(0); 3490 } 3491 3492 #undef __FUNCT__ 3493 #define __FUNCT__ "TSComputeRHSJacobianConstant" 3494 /*@C 3495 TSComputeRHSJacobianConstant - Reuses a Jacobian that is time-independent. 3496 3497 Collective on TS 3498 3499 Input Arguments: 3500 + ts - time stepping context 3501 . t - time at which to evaluate 3502 . U - state at which to evaluate 3503 - ctx - context 3504 3505 Output Arguments: 3506 + A - pointer to operator 3507 . B - pointer to preconditioning matrix 3508 - flg - matrix structure flag 3509 3510 Level: intermediate 3511 3512 Notes: 3513 This function is intended to be passed to TSSetRHSJacobian() to evaluate the Jacobian for linear time-independent problems. 3514 3515 .seealso: TSSetRHSFunction(), TSSetRHSJacobian(), TSComputeRHSFunctionLinear() 3516 @*/ 3517 PetscErrorCode TSComputeRHSJacobianConstant(TS ts,PetscReal t,Vec U,Mat *A,Mat *B,MatStructure *flg,void *ctx) 3518 { 3519 PetscFunctionBegin; 3520 *flg = SAME_PRECONDITIONER; 3521 PetscFunctionReturn(0); 3522 } 3523 3524 #undef __FUNCT__ 3525 #define __FUNCT__ "TSComputeIFunctionLinear" 3526 /*@C 3527 TSComputeIFunctionLinear - Evaluate the left hand side via the user-provided Jacobian, for linear problems only 3528 3529 Collective on TS 3530 3531 Input Arguments: 3532 + ts - time stepping context 3533 . t - time at which to evaluate 3534 . U - state at which to evaluate 3535 . Udot - time derivative of state vector 3536 - ctx - context 3537 3538 Output Arguments: 3539 . F - left hand side 3540 3541 Level: intermediate 3542 3543 Notes: 3544 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 3545 user is required to write their own TSComputeIFunction. 3546 This function is intended to be passed to TSSetIFunction() to evaluate the left hand side for linear problems. 3547 The matrix (and optionally the evaluation context) should be passed to TSSetIJacobian(). 3548 3549 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIJacobianConstant() 3550 @*/ 3551 PetscErrorCode TSComputeIFunctionLinear(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 3552 { 3553 PetscErrorCode ierr; 3554 Mat A,B; 3555 MatStructure flg2; 3556 3557 PetscFunctionBegin; 3558 ierr = TSGetIJacobian(ts,&A,&B,NULL,NULL);CHKERRQ(ierr); 3559 ierr = TSComputeIJacobian(ts,t,U,Udot,1.0,&A,&B,&flg2,PETSC_TRUE);CHKERRQ(ierr); 3560 ierr = MatMult(A,Udot,F);CHKERRQ(ierr); 3561 PetscFunctionReturn(0); 3562 } 3563 3564 #undef __FUNCT__ 3565 #define __FUNCT__ "TSComputeIJacobianConstant" 3566 /*@C 3567 TSComputeIJacobianConstant - Reuses a time-independent for a semi-implicit DAE or ODE 3568 3569 Collective on TS 3570 3571 Input Arguments: 3572 + ts - time stepping context 3573 . t - time at which to evaluate 3574 . U - state at which to evaluate 3575 . Udot - time derivative of state vector 3576 . shift - shift to apply 3577 - ctx - context 3578 3579 Output Arguments: 3580 + A - pointer to operator 3581 . B - pointer to preconditioning matrix 3582 - flg - matrix structure flag 3583 3584 Level: advanced 3585 3586 Notes: 3587 This function is intended to be passed to TSSetIJacobian() to evaluate the Jacobian for linear time-independent problems. 3588 3589 It is only appropriate for problems of the form 3590 3591 $ M Udot = F(U,t) 3592 3593 where M is constant and F is non-stiff. The user must pass M to TSSetIJacobian(). The current implementation only 3594 works with IMEX time integration methods such as TSROSW and TSARKIMEX, since there is no support for de-constructing 3595 an implicit operator of the form 3596 3597 $ shift*M + J 3598 3599 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 3600 a copy of M or reassemble it when requested. 3601 3602 .seealso: TSSetIFunction(), TSSetIJacobian(), TSComputeIFunctionLinear() 3603 @*/ 3604 PetscErrorCode TSComputeIJacobianConstant(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,void *ctx) 3605 { 3606 PetscErrorCode ierr; 3607 3608 PetscFunctionBegin; 3609 ierr = MatScale(*A, shift / ts->ijacobian.shift);CHKERRQ(ierr); 3610 ts->ijacobian.shift = shift; 3611 *flg = SAME_PRECONDITIONER; 3612 PetscFunctionReturn(0); 3613 } 3614 3615 #undef __FUNCT__ 3616 #define __FUNCT__ "TSGetEquationType" 3617 /*@ 3618 TSGetEquationType - Gets the type of the equation that TS is solving. 3619 3620 Not Collective 3621 3622 Input Parameter: 3623 . ts - the TS context 3624 3625 Output Parameter: 3626 . equation_type - see TSEquationType 3627 3628 Level: beginner 3629 3630 .keywords: TS, equation type 3631 3632 .seealso: TSSetEquationType(), TSEquationType 3633 @*/ 3634 PetscErrorCode TSGetEquationType(TS ts,TSEquationType *equation_type) 3635 { 3636 PetscFunctionBegin; 3637 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3638 PetscValidPointer(equation_type,2); 3639 *equation_type = ts->equation_type; 3640 PetscFunctionReturn(0); 3641 } 3642 3643 #undef __FUNCT__ 3644 #define __FUNCT__ "TSSetEquationType" 3645 /*@ 3646 TSSetEquationType - Sets the type of the equation that TS is solving. 3647 3648 Not Collective 3649 3650 Input Parameter: 3651 + ts - the TS context 3652 . equation_type - see TSEquationType 3653 3654 Level: advanced 3655 3656 .keywords: TS, equation type 3657 3658 .seealso: TSGetEquationType(), TSEquationType 3659 @*/ 3660 PetscErrorCode TSSetEquationType(TS ts,TSEquationType equation_type) 3661 { 3662 PetscFunctionBegin; 3663 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3664 ts->equation_type = equation_type; 3665 PetscFunctionReturn(0); 3666 } 3667 3668 #undef __FUNCT__ 3669 #define __FUNCT__ "TSGetConvergedReason" 3670 /*@ 3671 TSGetConvergedReason - Gets the reason the TS iteration was stopped. 3672 3673 Not Collective 3674 3675 Input Parameter: 3676 . ts - the TS context 3677 3678 Output Parameter: 3679 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3680 manual pages for the individual convergence tests for complete lists 3681 3682 Level: beginner 3683 3684 Notes: 3685 Can only be called after the call to TSSolve() is complete. 3686 3687 .keywords: TS, nonlinear, set, convergence, test 3688 3689 .seealso: TSSetConvergenceTest(), TSConvergedReason 3690 @*/ 3691 PetscErrorCode TSGetConvergedReason(TS ts,TSConvergedReason *reason) 3692 { 3693 PetscFunctionBegin; 3694 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3695 PetscValidPointer(reason,2); 3696 *reason = ts->reason; 3697 PetscFunctionReturn(0); 3698 } 3699 3700 #undef __FUNCT__ 3701 #define __FUNCT__ "TSSetConvergedReason" 3702 /*@ 3703 TSSetConvergedReason - Sets the reason for handling the convergence of TSSolve. 3704 3705 Not Collective 3706 3707 Input Parameter: 3708 + ts - the TS context 3709 . reason - negative value indicates diverged, positive value converged, see TSConvergedReason or the 3710 manual pages for the individual convergence tests for complete lists 3711 3712 Level: advanced 3713 3714 Notes: 3715 Can only be called during TSSolve() is active. 3716 3717 .keywords: TS, nonlinear, set, convergence, test 3718 3719 .seealso: TSConvergedReason 3720 @*/ 3721 PetscErrorCode TSSetConvergedReason(TS ts,TSConvergedReason reason) 3722 { 3723 PetscFunctionBegin; 3724 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3725 ts->reason = reason; 3726 PetscFunctionReturn(0); 3727 } 3728 3729 #undef __FUNCT__ 3730 #define __FUNCT__ "TSGetSolveTime" 3731 /*@ 3732 TSGetSolveTime - Gets the time after a call to TSSolve() 3733 3734 Not Collective 3735 3736 Input Parameter: 3737 . ts - the TS context 3738 3739 Output Parameter: 3740 . ftime - the final time. This time should correspond to the final time set with TSSetDuration() 3741 3742 Level: beginner 3743 3744 Notes: 3745 Can only be called after the call to TSSolve() is complete. 3746 3747 .keywords: TS, nonlinear, set, convergence, test 3748 3749 .seealso: TSSetConvergenceTest(), TSConvergedReason 3750 @*/ 3751 PetscErrorCode TSGetSolveTime(TS ts,PetscReal *ftime) 3752 { 3753 PetscFunctionBegin; 3754 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3755 PetscValidPointer(ftime,2); 3756 *ftime = ts->solvetime; 3757 PetscFunctionReturn(0); 3758 } 3759 3760 #undef __FUNCT__ 3761 #define __FUNCT__ "TSGetSNESIterations" 3762 /*@ 3763 TSGetSNESIterations - Gets the total number of nonlinear iterations 3764 used by the time integrator. 3765 3766 Not Collective 3767 3768 Input Parameter: 3769 . ts - TS context 3770 3771 Output Parameter: 3772 . nits - number of nonlinear iterations 3773 3774 Notes: 3775 This counter is reset to zero for each successive call to TSSolve(). 3776 3777 Level: intermediate 3778 3779 .keywords: TS, get, number, nonlinear, iterations 3780 3781 .seealso: TSGetKSPIterations() 3782 @*/ 3783 PetscErrorCode TSGetSNESIterations(TS ts,PetscInt *nits) 3784 { 3785 PetscFunctionBegin; 3786 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3787 PetscValidIntPointer(nits,2); 3788 *nits = ts->snes_its; 3789 PetscFunctionReturn(0); 3790 } 3791 3792 #undef __FUNCT__ 3793 #define __FUNCT__ "TSGetKSPIterations" 3794 /*@ 3795 TSGetKSPIterations - Gets the total number of linear iterations 3796 used by the time integrator. 3797 3798 Not Collective 3799 3800 Input Parameter: 3801 . ts - TS context 3802 3803 Output Parameter: 3804 . lits - number of linear iterations 3805 3806 Notes: 3807 This counter is reset to zero for each successive call to TSSolve(). 3808 3809 Level: intermediate 3810 3811 .keywords: TS, get, number, linear, iterations 3812 3813 .seealso: TSGetSNESIterations(), SNESGetKSPIterations() 3814 @*/ 3815 PetscErrorCode TSGetKSPIterations(TS ts,PetscInt *lits) 3816 { 3817 PetscFunctionBegin; 3818 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3819 PetscValidIntPointer(lits,2); 3820 *lits = ts->ksp_its; 3821 PetscFunctionReturn(0); 3822 } 3823 3824 #undef __FUNCT__ 3825 #define __FUNCT__ "TSGetStepRejections" 3826 /*@ 3827 TSGetStepRejections - Gets the total number of rejected steps. 3828 3829 Not Collective 3830 3831 Input Parameter: 3832 . ts - TS context 3833 3834 Output Parameter: 3835 . rejects - number of steps rejected 3836 3837 Notes: 3838 This counter is reset to zero for each successive call to TSSolve(). 3839 3840 Level: intermediate 3841 3842 .keywords: TS, get, number 3843 3844 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetSNESFailures(), TSSetMaxSNESFailures(), TSSetErrorIfStepFails() 3845 @*/ 3846 PetscErrorCode TSGetStepRejections(TS ts,PetscInt *rejects) 3847 { 3848 PetscFunctionBegin; 3849 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3850 PetscValidIntPointer(rejects,2); 3851 *rejects = ts->reject; 3852 PetscFunctionReturn(0); 3853 } 3854 3855 #undef __FUNCT__ 3856 #define __FUNCT__ "TSGetSNESFailures" 3857 /*@ 3858 TSGetSNESFailures - Gets the total number of failed SNES solves 3859 3860 Not Collective 3861 3862 Input Parameter: 3863 . ts - TS context 3864 3865 Output Parameter: 3866 . fails - number of failed nonlinear solves 3867 3868 Notes: 3869 This counter is reset to zero for each successive call to TSSolve(). 3870 3871 Level: intermediate 3872 3873 .keywords: TS, get, number 3874 3875 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSSetMaxSNESFailures() 3876 @*/ 3877 PetscErrorCode TSGetSNESFailures(TS ts,PetscInt *fails) 3878 { 3879 PetscFunctionBegin; 3880 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3881 PetscValidIntPointer(fails,2); 3882 *fails = ts->num_snes_failures; 3883 PetscFunctionReturn(0); 3884 } 3885 3886 #undef __FUNCT__ 3887 #define __FUNCT__ "TSSetMaxStepRejections" 3888 /*@ 3889 TSSetMaxStepRejections - Sets the maximum number of step rejections before a step fails 3890 3891 Not Collective 3892 3893 Input Parameter: 3894 + ts - TS context 3895 - rejects - maximum number of rejected steps, pass -1 for unlimited 3896 3897 Notes: 3898 The counter is reset to zero for each step 3899 3900 Options Database Key: 3901 . -ts_max_reject - Maximum number of step rejections before a step fails 3902 3903 Level: intermediate 3904 3905 .keywords: TS, set, maximum, number 3906 3907 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxSNESFailures(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 3908 @*/ 3909 PetscErrorCode TSSetMaxStepRejections(TS ts,PetscInt rejects) 3910 { 3911 PetscFunctionBegin; 3912 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3913 ts->max_reject = rejects; 3914 PetscFunctionReturn(0); 3915 } 3916 3917 #undef __FUNCT__ 3918 #define __FUNCT__ "TSSetMaxSNESFailures" 3919 /*@ 3920 TSSetMaxSNESFailures - Sets the maximum number of failed SNES solves 3921 3922 Not Collective 3923 3924 Input Parameter: 3925 + ts - TS context 3926 - fails - maximum number of failed nonlinear solves, pass -1 for unlimited 3927 3928 Notes: 3929 The counter is reset to zero for each successive call to TSSolve(). 3930 3931 Options Database Key: 3932 . -ts_max_snes_failures - Maximum number of nonlinear solve failures 3933 3934 Level: intermediate 3935 3936 .keywords: TS, set, maximum, number 3937 3938 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), SNESGetConvergedReason(), TSGetConvergedReason() 3939 @*/ 3940 PetscErrorCode TSSetMaxSNESFailures(TS ts,PetscInt fails) 3941 { 3942 PetscFunctionBegin; 3943 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3944 ts->max_snes_failures = fails; 3945 PetscFunctionReturn(0); 3946 } 3947 3948 #undef __FUNCT__ 3949 #define __FUNCT__ "TSSetErrorIfStepFails()" 3950 /*@ 3951 TSSetErrorIfStepFails - Error if no step succeeds 3952 3953 Not Collective 3954 3955 Input Parameter: 3956 + ts - TS context 3957 - err - PETSC_TRUE to error if no step succeeds, PETSC_FALSE to return without failure 3958 3959 Options Database Key: 3960 . -ts_error_if_step_fails - Error if no step succeeds 3961 3962 Level: intermediate 3963 3964 .keywords: TS, set, error 3965 3966 .seealso: TSGetSNESIterations(), TSGetKSPIterations(), TSSetMaxStepRejections(), TSGetStepRejections(), TSGetSNESFailures(), TSSetErrorIfStepFails(), TSGetConvergedReason() 3967 @*/ 3968 PetscErrorCode TSSetErrorIfStepFails(TS ts,PetscBool err) 3969 { 3970 PetscFunctionBegin; 3971 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 3972 ts->errorifstepfailed = err; 3973 PetscFunctionReturn(0); 3974 } 3975 3976 #undef __FUNCT__ 3977 #define __FUNCT__ "TSMonitorSolutionBinary" 3978 /*@C 3979 TSMonitorSolutionBinary - Monitors progress of the TS solvers by VecView() for the solution at each timestep. Normally the viewer is a binary file 3980 3981 Collective on TS 3982 3983 Input Parameters: 3984 + ts - the TS context 3985 . step - current time-step 3986 . ptime - current time 3987 . u - current state 3988 - viewer - binary viewer 3989 3990 Level: intermediate 3991 3992 .keywords: TS, vector, monitor, view 3993 3994 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 3995 @*/ 3996 PetscErrorCode TSMonitorSolutionBinary(TS ts,PetscInt step,PetscReal ptime,Vec u,void *viewer) 3997 { 3998 PetscErrorCode ierr; 3999 PetscViewer v = (PetscViewer)viewer; 4000 4001 PetscFunctionBegin; 4002 ierr = VecView(u,v);CHKERRQ(ierr); 4003 PetscFunctionReturn(0); 4004 } 4005 4006 #undef __FUNCT__ 4007 #define __FUNCT__ "TSMonitorSolutionVTK" 4008 /*@C 4009 TSMonitorSolutionVTK - Monitors progress of the TS solvers by VecView() for the solution at each timestep. 4010 4011 Collective on TS 4012 4013 Input Parameters: 4014 + ts - the TS context 4015 . step - current time-step 4016 . ptime - current time 4017 . u - current state 4018 - filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4019 4020 Level: intermediate 4021 4022 Notes: 4023 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. 4024 These are named according to the file name template. 4025 4026 This function is normally passed as an argument to TSMonitorSet() along with TSMonitorSolutionVTKDestroy(). 4027 4028 .keywords: TS, vector, monitor, view 4029 4030 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4031 @*/ 4032 PetscErrorCode TSMonitorSolutionVTK(TS ts,PetscInt step,PetscReal ptime,Vec u,void *filenametemplate) 4033 { 4034 PetscErrorCode ierr; 4035 char filename[PETSC_MAX_PATH_LEN]; 4036 PetscViewer viewer; 4037 4038 PetscFunctionBegin; 4039 ierr = PetscSNPrintf(filename,sizeof(filename),(const char*)filenametemplate,step);CHKERRQ(ierr); 4040 ierr = PetscViewerVTKOpen(PetscObjectComm((PetscObject)ts),filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 4041 ierr = VecView(u,viewer);CHKERRQ(ierr); 4042 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 4043 PetscFunctionReturn(0); 4044 } 4045 4046 #undef __FUNCT__ 4047 #define __FUNCT__ "TSMonitorSolutionVTKDestroy" 4048 /*@C 4049 TSMonitorSolutionVTKDestroy - Destroy context for monitoring 4050 4051 Collective on TS 4052 4053 Input Parameters: 4054 . filenametemplate - string containing a format specifier for the integer time step (e.g. %03D) 4055 4056 Level: intermediate 4057 4058 Note: 4059 This function is normally passed to TSMonitorSet() along with TSMonitorSolutionVTK(). 4060 4061 .keywords: TS, vector, monitor, view 4062 4063 .seealso: TSMonitorSet(), TSMonitorSolutionVTK() 4064 @*/ 4065 PetscErrorCode TSMonitorSolutionVTKDestroy(void *filenametemplate) 4066 { 4067 PetscErrorCode ierr; 4068 4069 PetscFunctionBegin; 4070 ierr = PetscFree(*(char**)filenametemplate);CHKERRQ(ierr); 4071 PetscFunctionReturn(0); 4072 } 4073 4074 #undef __FUNCT__ 4075 #define __FUNCT__ "TSGetAdapt" 4076 /*@ 4077 TSGetAdapt - Get the adaptive controller context for the current method 4078 4079 Collective on TS if controller has not been created yet 4080 4081 Input Arguments: 4082 . ts - time stepping context 4083 4084 Output Arguments: 4085 . adapt - adaptive controller 4086 4087 Level: intermediate 4088 4089 .seealso: TSAdapt, TSAdaptSetType(), TSAdaptChoose() 4090 @*/ 4091 PetscErrorCode TSGetAdapt(TS ts,TSAdapt *adapt) 4092 { 4093 PetscErrorCode ierr; 4094 4095 PetscFunctionBegin; 4096 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4097 PetscValidPointer(adapt,2); 4098 if (!ts->adapt) { 4099 ierr = TSAdaptCreate(PetscObjectComm((PetscObject)ts),&ts->adapt);CHKERRQ(ierr); 4100 ierr = PetscLogObjectParent(ts,ts->adapt);CHKERRQ(ierr); 4101 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->adapt,(PetscObject)ts,1);CHKERRQ(ierr); 4102 } 4103 *adapt = ts->adapt; 4104 PetscFunctionReturn(0); 4105 } 4106 4107 #undef __FUNCT__ 4108 #define __FUNCT__ "TSSetTolerances" 4109 /*@ 4110 TSSetTolerances - Set tolerances for local truncation error when using adaptive controller 4111 4112 Logically Collective 4113 4114 Input Arguments: 4115 + ts - time integration context 4116 . atol - scalar absolute tolerances, PETSC_DECIDE to leave current value 4117 . vatol - vector of absolute tolerances or NULL, used in preference to atol if present 4118 . rtol - scalar relative tolerances, PETSC_DECIDE to leave current value 4119 - vrtol - vector of relative tolerances or NULL, used in preference to atol if present 4120 4121 Level: beginner 4122 4123 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSGetTolerances() 4124 @*/ 4125 PetscErrorCode TSSetTolerances(TS ts,PetscReal atol,Vec vatol,PetscReal rtol,Vec vrtol) 4126 { 4127 PetscErrorCode ierr; 4128 4129 PetscFunctionBegin; 4130 if (atol != PETSC_DECIDE && atol != PETSC_DEFAULT) ts->atol = atol; 4131 if (vatol) { 4132 ierr = PetscObjectReference((PetscObject)vatol);CHKERRQ(ierr); 4133 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 4134 4135 ts->vatol = vatol; 4136 } 4137 if (rtol != PETSC_DECIDE && rtol != PETSC_DEFAULT) ts->rtol = rtol; 4138 if (vrtol) { 4139 ierr = PetscObjectReference((PetscObject)vrtol);CHKERRQ(ierr); 4140 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 4141 4142 ts->vrtol = vrtol; 4143 } 4144 PetscFunctionReturn(0); 4145 } 4146 4147 #undef __FUNCT__ 4148 #define __FUNCT__ "TSGetTolerances" 4149 /*@ 4150 TSGetTolerances - Get tolerances for local truncation error when using adaptive controller 4151 4152 Logically Collective 4153 4154 Input Arguments: 4155 . ts - time integration context 4156 4157 Output Arguments: 4158 + atol - scalar absolute tolerances, NULL to ignore 4159 . vatol - vector of absolute tolerances, NULL to ignore 4160 . rtol - scalar relative tolerances, NULL to ignore 4161 - vrtol - vector of relative tolerances, NULL to ignore 4162 4163 Level: beginner 4164 4165 .seealso: TS, TSAdapt, TSVecNormWRMS(), TSSetTolerances() 4166 @*/ 4167 PetscErrorCode TSGetTolerances(TS ts,PetscReal *atol,Vec *vatol,PetscReal *rtol,Vec *vrtol) 4168 { 4169 PetscFunctionBegin; 4170 if (atol) *atol = ts->atol; 4171 if (vatol) *vatol = ts->vatol; 4172 if (rtol) *rtol = ts->rtol; 4173 if (vrtol) *vrtol = ts->vrtol; 4174 PetscFunctionReturn(0); 4175 } 4176 4177 #undef __FUNCT__ 4178 #define __FUNCT__ "TSErrorNormWRMS" 4179 /*@ 4180 TSErrorNormWRMS - compute a weighted norm of the difference between a vector and the current state 4181 4182 Collective on TS 4183 4184 Input Arguments: 4185 + ts - time stepping context 4186 - Y - state vector to be compared to ts->vec_sol 4187 4188 Output Arguments: 4189 . norm - weighted norm, a value of 1.0 is considered small 4190 4191 Level: developer 4192 4193 .seealso: TSSetTolerances() 4194 @*/ 4195 PetscErrorCode TSErrorNormWRMS(TS ts,Vec Y,PetscReal *norm) 4196 { 4197 PetscErrorCode ierr; 4198 PetscInt i,n,N; 4199 const PetscScalar *u,*y; 4200 Vec U; 4201 PetscReal sum,gsum; 4202 4203 PetscFunctionBegin; 4204 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4205 PetscValidHeaderSpecific(Y,VEC_CLASSID,2); 4206 PetscValidPointer(norm,3); 4207 U = ts->vec_sol; 4208 PetscCheckSameTypeAndComm(U,1,Y,2); 4209 if (U == Y) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_IDN,"Y cannot be the TS solution vector"); 4210 4211 ierr = VecGetSize(U,&N);CHKERRQ(ierr); 4212 ierr = VecGetLocalSize(U,&n);CHKERRQ(ierr); 4213 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 4214 ierr = VecGetArrayRead(Y,&y);CHKERRQ(ierr); 4215 sum = 0.; 4216 if (ts->vatol && ts->vrtol) { 4217 const PetscScalar *atol,*rtol; 4218 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4219 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4220 for (i=0; i<n; i++) { 4221 PetscReal tol = PetscRealPart(atol[i]) + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4222 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4223 } 4224 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4225 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4226 } else if (ts->vatol) { /* vector atol, scalar rtol */ 4227 const PetscScalar *atol; 4228 ierr = VecGetArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4229 for (i=0; i<n; i++) { 4230 PetscReal tol = PetscRealPart(atol[i]) + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4231 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4232 } 4233 ierr = VecRestoreArrayRead(ts->vatol,&atol);CHKERRQ(ierr); 4234 } else if (ts->vrtol) { /* scalar atol, vector rtol */ 4235 const PetscScalar *rtol; 4236 ierr = VecGetArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4237 for (i=0; i<n; i++) { 4238 PetscReal tol = ts->atol + PetscRealPart(rtol[i]) * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4239 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4240 } 4241 ierr = VecRestoreArrayRead(ts->vrtol,&rtol);CHKERRQ(ierr); 4242 } else { /* scalar atol, scalar rtol */ 4243 for (i=0; i<n; i++) { 4244 PetscReal tol = ts->atol + ts->rtol * PetscMax(PetscAbsScalar(u[i]),PetscAbsScalar(y[i])); 4245 sum += PetscSqr(PetscAbsScalar(y[i] - u[i]) / tol); 4246 } 4247 } 4248 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 4249 ierr = VecRestoreArrayRead(Y,&y);CHKERRQ(ierr); 4250 4251 ierr = MPI_Allreduce(&sum,&gsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4252 *norm = PetscSqrtReal(gsum / N); 4253 if (PetscIsInfOrNanScalar(*norm)) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 4254 PetscFunctionReturn(0); 4255 } 4256 4257 #undef __FUNCT__ 4258 #define __FUNCT__ "TSSetCFLTimeLocal" 4259 /*@ 4260 TSSetCFLTimeLocal - Set the local CFL constraint relative to forward Euler 4261 4262 Logically Collective on TS 4263 4264 Input Arguments: 4265 + ts - time stepping context 4266 - cfltime - maximum stable time step if using forward Euler (value can be different on each process) 4267 4268 Note: 4269 After calling this function, the global CFL time can be obtained by calling TSGetCFLTime() 4270 4271 Level: intermediate 4272 4273 .seealso: TSGetCFLTime(), TSADAPTCFL 4274 @*/ 4275 PetscErrorCode TSSetCFLTimeLocal(TS ts,PetscReal cfltime) 4276 { 4277 PetscFunctionBegin; 4278 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4279 ts->cfltime_local = cfltime; 4280 ts->cfltime = -1.; 4281 PetscFunctionReturn(0); 4282 } 4283 4284 #undef __FUNCT__ 4285 #define __FUNCT__ "TSGetCFLTime" 4286 /*@ 4287 TSGetCFLTime - Get the maximum stable time step according to CFL criteria applied to forward Euler 4288 4289 Collective on TS 4290 4291 Input Arguments: 4292 . ts - time stepping context 4293 4294 Output Arguments: 4295 . cfltime - maximum stable time step for forward Euler 4296 4297 Level: advanced 4298 4299 .seealso: TSSetCFLTimeLocal() 4300 @*/ 4301 PetscErrorCode TSGetCFLTime(TS ts,PetscReal *cfltime) 4302 { 4303 PetscErrorCode ierr; 4304 4305 PetscFunctionBegin; 4306 if (ts->cfltime < 0) { 4307 ierr = MPI_Allreduce(&ts->cfltime_local,&ts->cfltime,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)ts));CHKERRQ(ierr); 4308 } 4309 *cfltime = ts->cfltime; 4310 PetscFunctionReturn(0); 4311 } 4312 4313 #undef __FUNCT__ 4314 #define __FUNCT__ "TSVISetVariableBounds" 4315 /*@ 4316 TSVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu 4317 4318 Input Parameters: 4319 . ts - the TS context. 4320 . xl - lower bound. 4321 . xu - upper bound. 4322 4323 Notes: 4324 If this routine is not called then the lower and upper bounds are set to 4325 SNES_VI_NINF and SNES_VI_INF respectively during SNESSetUp(). 4326 4327 Level: advanced 4328 4329 @*/ 4330 PetscErrorCode TSVISetVariableBounds(TS ts, Vec xl, Vec xu) 4331 { 4332 PetscErrorCode ierr; 4333 SNES snes; 4334 4335 PetscFunctionBegin; 4336 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 4337 ierr = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(ierr); 4338 PetscFunctionReturn(0); 4339 } 4340 4341 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4342 #include <mex.h> 4343 4344 typedef struct {char *funcname; mxArray *ctx;} TSMatlabContext; 4345 4346 #undef __FUNCT__ 4347 #define __FUNCT__ "TSComputeFunction_Matlab" 4348 /* 4349 TSComputeFunction_Matlab - Calls the function that has been set with 4350 TSSetFunctionMatlab(). 4351 4352 Collective on TS 4353 4354 Input Parameters: 4355 + snes - the TS context 4356 - u - input vector 4357 4358 Output Parameter: 4359 . y - function vector, as set by TSSetFunction() 4360 4361 Notes: 4362 TSComputeFunction() is typically used within nonlinear solvers 4363 implementations, so most users would not generally call this routine 4364 themselves. 4365 4366 Level: developer 4367 4368 .keywords: TS, nonlinear, compute, function 4369 4370 .seealso: TSSetFunction(), TSGetFunction() 4371 */ 4372 PetscErrorCode TSComputeFunction_Matlab(TS snes,PetscReal time,Vec u,Vec udot,Vec y, void *ctx) 4373 { 4374 PetscErrorCode ierr; 4375 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4376 int nlhs = 1,nrhs = 7; 4377 mxArray *plhs[1],*prhs[7]; 4378 long long int lx = 0,lxdot = 0,ly = 0,ls = 0; 4379 4380 PetscFunctionBegin; 4381 PetscValidHeaderSpecific(snes,TS_CLASSID,1); 4382 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4383 PetscValidHeaderSpecific(udot,VEC_CLASSID,4); 4384 PetscValidHeaderSpecific(y,VEC_CLASSID,5); 4385 PetscCheckSameComm(snes,1,u,3); 4386 PetscCheckSameComm(snes,1,y,5); 4387 4388 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4389 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4390 ierr = PetscMemcpy(&lxdot,&udot,sizeof(udot));CHKERRQ(ierr); 4391 ierr = PetscMemcpy(&ly,&y,sizeof(u));CHKERRQ(ierr); 4392 4393 prhs[0] = mxCreateDoubleScalar((double)ls); 4394 prhs[1] = mxCreateDoubleScalar(time); 4395 prhs[2] = mxCreateDoubleScalar((double)lx); 4396 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4397 prhs[4] = mxCreateDoubleScalar((double)ly); 4398 prhs[5] = mxCreateString(sctx->funcname); 4399 prhs[6] = sctx->ctx; 4400 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeFunctionInternal");CHKERRQ(ierr); 4401 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4402 mxDestroyArray(prhs[0]); 4403 mxDestroyArray(prhs[1]); 4404 mxDestroyArray(prhs[2]); 4405 mxDestroyArray(prhs[3]); 4406 mxDestroyArray(prhs[4]); 4407 mxDestroyArray(prhs[5]); 4408 mxDestroyArray(plhs[0]); 4409 PetscFunctionReturn(0); 4410 } 4411 4412 4413 #undef __FUNCT__ 4414 #define __FUNCT__ "TSSetFunctionMatlab" 4415 /* 4416 TSSetFunctionMatlab - Sets the function evaluation routine and function 4417 vector for use by the TS routines in solving ODEs 4418 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4419 4420 Logically Collective on TS 4421 4422 Input Parameters: 4423 + ts - the TS context 4424 - func - function evaluation routine 4425 4426 Calling sequence of func: 4427 $ func (TS ts,PetscReal time,Vec u,Vec udot,Vec f,void *ctx); 4428 4429 Level: beginner 4430 4431 .keywords: TS, nonlinear, set, function 4432 4433 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4434 */ 4435 PetscErrorCode TSSetFunctionMatlab(TS ts,const char *func,mxArray *ctx) 4436 { 4437 PetscErrorCode ierr; 4438 TSMatlabContext *sctx; 4439 4440 PetscFunctionBegin; 4441 /* currently sctx is memory bleed */ 4442 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4443 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4444 /* 4445 This should work, but it doesn't 4446 sctx->ctx = ctx; 4447 mexMakeArrayPersistent(sctx->ctx); 4448 */ 4449 sctx->ctx = mxDuplicateArray(ctx); 4450 4451 ierr = TSSetIFunction(ts,NULL,TSComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4452 PetscFunctionReturn(0); 4453 } 4454 4455 #undef __FUNCT__ 4456 #define __FUNCT__ "TSComputeJacobian_Matlab" 4457 /* 4458 TSComputeJacobian_Matlab - Calls the function that has been set with 4459 TSSetJacobianMatlab(). 4460 4461 Collective on TS 4462 4463 Input Parameters: 4464 + ts - the TS context 4465 . u - input vector 4466 . A, B - the matrices 4467 - ctx - user context 4468 4469 Output Parameter: 4470 . flag - structure of the matrix 4471 4472 Level: developer 4473 4474 .keywords: TS, nonlinear, compute, function 4475 4476 .seealso: TSSetFunction(), TSGetFunction() 4477 @*/ 4478 PetscErrorCode TSComputeJacobian_Matlab(TS ts,PetscReal time,Vec u,Vec udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4479 { 4480 PetscErrorCode ierr; 4481 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4482 int nlhs = 2,nrhs = 9; 4483 mxArray *plhs[2],*prhs[9]; 4484 long long int lx = 0,lxdot = 0,lA = 0,ls = 0, lB = 0; 4485 4486 PetscFunctionBegin; 4487 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4488 PetscValidHeaderSpecific(u,VEC_CLASSID,3); 4489 4490 /* call Matlab function in ctx with arguments u and y */ 4491 4492 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4493 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4494 ierr = PetscMemcpy(&lxdot,&udot,sizeof(u));CHKERRQ(ierr); 4495 ierr = PetscMemcpy(&lA,A,sizeof(u));CHKERRQ(ierr); 4496 ierr = PetscMemcpy(&lB,B,sizeof(u));CHKERRQ(ierr); 4497 4498 prhs[0] = mxCreateDoubleScalar((double)ls); 4499 prhs[1] = mxCreateDoubleScalar((double)time); 4500 prhs[2] = mxCreateDoubleScalar((double)lx); 4501 prhs[3] = mxCreateDoubleScalar((double)lxdot); 4502 prhs[4] = mxCreateDoubleScalar((double)shift); 4503 prhs[5] = mxCreateDoubleScalar((double)lA); 4504 prhs[6] = mxCreateDoubleScalar((double)lB); 4505 prhs[7] = mxCreateString(sctx->funcname); 4506 prhs[8] = sctx->ctx; 4507 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSComputeJacobianInternal");CHKERRQ(ierr); 4508 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4509 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4510 mxDestroyArray(prhs[0]); 4511 mxDestroyArray(prhs[1]); 4512 mxDestroyArray(prhs[2]); 4513 mxDestroyArray(prhs[3]); 4514 mxDestroyArray(prhs[4]); 4515 mxDestroyArray(prhs[5]); 4516 mxDestroyArray(prhs[6]); 4517 mxDestroyArray(prhs[7]); 4518 mxDestroyArray(plhs[0]); 4519 mxDestroyArray(plhs[1]); 4520 PetscFunctionReturn(0); 4521 } 4522 4523 4524 #undef __FUNCT__ 4525 #define __FUNCT__ "TSSetJacobianMatlab" 4526 /* 4527 TSSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4528 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 4529 4530 Logically Collective on TS 4531 4532 Input Parameters: 4533 + ts - the TS context 4534 . A,B - Jacobian matrices 4535 . func - function evaluation routine 4536 - ctx - user context 4537 4538 Calling sequence of func: 4539 $ flag = func (TS ts,PetscReal time,Vec u,Vec udot,Mat A,Mat B,void *ctx); 4540 4541 4542 Level: developer 4543 4544 .keywords: TS, nonlinear, set, function 4545 4546 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4547 */ 4548 PetscErrorCode TSSetJacobianMatlab(TS ts,Mat A,Mat B,const char *func,mxArray *ctx) 4549 { 4550 PetscErrorCode ierr; 4551 TSMatlabContext *sctx; 4552 4553 PetscFunctionBegin; 4554 /* currently sctx is memory bleed */ 4555 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4556 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4557 /* 4558 This should work, but it doesn't 4559 sctx->ctx = ctx; 4560 mexMakeArrayPersistent(sctx->ctx); 4561 */ 4562 sctx->ctx = mxDuplicateArray(ctx); 4563 4564 ierr = TSSetIJacobian(ts,A,B,TSComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4565 PetscFunctionReturn(0); 4566 } 4567 4568 #undef __FUNCT__ 4569 #define __FUNCT__ "TSMonitor_Matlab" 4570 /* 4571 TSMonitor_Matlab - Calls the function that has been set with TSMonitorSetMatlab(). 4572 4573 Collective on TS 4574 4575 .seealso: TSSetFunction(), TSGetFunction() 4576 @*/ 4577 PetscErrorCode TSMonitor_Matlab(TS ts,PetscInt it, PetscReal time,Vec u, void *ctx) 4578 { 4579 PetscErrorCode ierr; 4580 TSMatlabContext *sctx = (TSMatlabContext*)ctx; 4581 int nlhs = 1,nrhs = 6; 4582 mxArray *plhs[1],*prhs[6]; 4583 long long int lx = 0,ls = 0; 4584 4585 PetscFunctionBegin; 4586 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4587 PetscValidHeaderSpecific(u,VEC_CLASSID,4); 4588 4589 ierr = PetscMemcpy(&ls,&ts,sizeof(ts));CHKERRQ(ierr); 4590 ierr = PetscMemcpy(&lx,&u,sizeof(u));CHKERRQ(ierr); 4591 4592 prhs[0] = mxCreateDoubleScalar((double)ls); 4593 prhs[1] = mxCreateDoubleScalar((double)it); 4594 prhs[2] = mxCreateDoubleScalar((double)time); 4595 prhs[3] = mxCreateDoubleScalar((double)lx); 4596 prhs[4] = mxCreateString(sctx->funcname); 4597 prhs[5] = sctx->ctx; 4598 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscTSMonitorInternal");CHKERRQ(ierr); 4599 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4600 mxDestroyArray(prhs[0]); 4601 mxDestroyArray(prhs[1]); 4602 mxDestroyArray(prhs[2]); 4603 mxDestroyArray(prhs[3]); 4604 mxDestroyArray(prhs[4]); 4605 mxDestroyArray(plhs[0]); 4606 PetscFunctionReturn(0); 4607 } 4608 4609 4610 #undef __FUNCT__ 4611 #define __FUNCT__ "TSMonitorSetMatlab" 4612 /* 4613 TSMonitorSetMatlab - Sets the monitor function from Matlab 4614 4615 Level: developer 4616 4617 .keywords: TS, nonlinear, set, function 4618 4619 .seealso: TSGetFunction(), TSComputeFunction(), TSSetJacobian(), TSSetFunction() 4620 */ 4621 PetscErrorCode TSMonitorSetMatlab(TS ts,const char *func,mxArray *ctx) 4622 { 4623 PetscErrorCode ierr; 4624 TSMatlabContext *sctx; 4625 4626 PetscFunctionBegin; 4627 /* currently sctx is memory bleed */ 4628 ierr = PetscMalloc(sizeof(TSMatlabContext),&sctx);CHKERRQ(ierr); 4629 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4630 /* 4631 This should work, but it doesn't 4632 sctx->ctx = ctx; 4633 mexMakeArrayPersistent(sctx->ctx); 4634 */ 4635 sctx->ctx = mxDuplicateArray(ctx); 4636 4637 ierr = TSMonitorSet(ts,TSMonitor_Matlab,sctx,NULL);CHKERRQ(ierr); 4638 PetscFunctionReturn(0); 4639 } 4640 #endif 4641 4642 4643 4644 #undef __FUNCT__ 4645 #define __FUNCT__ "TSMonitorLGSolution" 4646 /*@C 4647 TSMonitorLGSolution - Monitors progress of the TS solvers by plotting each component of the solution vector 4648 in a time based line graph 4649 4650 Collective on TS 4651 4652 Input Parameters: 4653 + ts - the TS context 4654 . step - current time-step 4655 . ptime - current time 4656 - lg - a line graph object 4657 4658 Level: intermediate 4659 4660 Notes: each process in a parallel run displays its component solutions in a separate window 4661 4662 .keywords: TS, vector, monitor, view 4663 4664 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView() 4665 @*/ 4666 PetscErrorCode TSMonitorLGSolution(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4667 { 4668 PetscErrorCode ierr; 4669 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4670 const PetscScalar *yy; 4671 PetscInt dim; 4672 4673 PetscFunctionBegin; 4674 if (!step) { 4675 PetscDrawAxis axis; 4676 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4677 ierr = PetscDrawAxisSetLabels(axis,"Solution as function of time","Time","Solution");CHKERRQ(ierr); 4678 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4679 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4680 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4681 } 4682 ierr = VecGetArrayRead(u,&yy);CHKERRQ(ierr); 4683 #if defined(PETSC_USE_COMPLEX) 4684 { 4685 PetscReal *yreal; 4686 PetscInt i,n; 4687 ierr = VecGetLocalSize(u,&n);CHKERRQ(ierr); 4688 ierr = PetscMalloc(n*sizeof(PetscReal),&yreal);CHKERRQ(ierr); 4689 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4690 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4691 ierr = PetscFree(yreal);CHKERRQ(ierr); 4692 } 4693 #else 4694 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4695 #endif 4696 ierr = VecRestoreArrayRead(u,&yy);CHKERRQ(ierr); 4697 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 4698 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4699 } 4700 PetscFunctionReturn(0); 4701 } 4702 4703 #undef __FUNCT__ 4704 #define __FUNCT__ "TSMonitorLGError" 4705 /*@C 4706 TSMonitorLGError - Monitors progress of the TS solvers by plotting each component of the solution vector 4707 in a time based line graph 4708 4709 Collective on TS 4710 4711 Input Parameters: 4712 + ts - the TS context 4713 . step - current time-step 4714 . ptime - current time 4715 - lg - a line graph object 4716 4717 Level: intermediate 4718 4719 Notes: 4720 Only for sequential solves. 4721 4722 The user must provide the solution using TSSetSolutionFunction() to use this monitor. 4723 4724 Options Database Keys: 4725 . -ts_monitor_lg_error - create a graphical monitor of error history 4726 4727 .keywords: TS, vector, monitor, view 4728 4729 .seealso: TSMonitorSet(), TSMonitorDefault(), VecView(), TSSetSolutionFunction() 4730 @*/ 4731 PetscErrorCode TSMonitorLGError(TS ts,PetscInt step,PetscReal ptime,Vec u,void *dummy) 4732 { 4733 PetscErrorCode ierr; 4734 TSMonitorLGCtx ctx = (TSMonitorLGCtx)dummy; 4735 const PetscScalar *yy; 4736 Vec y; 4737 PetscInt dim; 4738 4739 PetscFunctionBegin; 4740 if (!step) { 4741 PetscDrawAxis axis; 4742 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4743 ierr = PetscDrawAxisSetLabels(axis,"Error in solution as function of time","Time","Solution");CHKERRQ(ierr); 4744 ierr = VecGetLocalSize(u,&dim);CHKERRQ(ierr); 4745 ierr = PetscDrawLGSetDimension(ctx->lg,dim);CHKERRQ(ierr); 4746 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4747 } 4748 ierr = VecDuplicate(u,&y);CHKERRQ(ierr); 4749 ierr = TSComputeSolutionFunction(ts,ptime,y);CHKERRQ(ierr); 4750 ierr = VecAXPY(y,-1.0,u);CHKERRQ(ierr); 4751 ierr = VecGetArrayRead(y,&yy);CHKERRQ(ierr); 4752 #if defined(PETSC_USE_COMPLEX) 4753 { 4754 PetscReal *yreal; 4755 PetscInt i,n; 4756 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 4757 ierr = PetscMalloc(n*sizeof(PetscReal),&yreal);CHKERRQ(ierr); 4758 for (i=0; i<n; i++) yreal[i] = PetscRealPart(yy[i]); 4759 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yreal);CHKERRQ(ierr); 4760 ierr = PetscFree(yreal);CHKERRQ(ierr); 4761 } 4762 #else 4763 ierr = PetscDrawLGAddCommonPoint(ctx->lg,ptime,yy);CHKERRQ(ierr); 4764 #endif 4765 ierr = VecRestoreArrayRead(y,&yy);CHKERRQ(ierr); 4766 ierr = VecDestroy(&y);CHKERRQ(ierr); 4767 if (((ctx->howoften > 0) && (!(step % ctx->howoften))) || ((ctx->howoften == -1) && ts->reason)) { 4768 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4769 } 4770 PetscFunctionReturn(0); 4771 } 4772 4773 #undef __FUNCT__ 4774 #define __FUNCT__ "TSMonitorLGSNESIterations" 4775 PetscErrorCode TSMonitorLGSNESIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4776 { 4777 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4778 PetscReal x = ptime,y; 4779 PetscErrorCode ierr; 4780 PetscInt its; 4781 4782 PetscFunctionBegin; 4783 if (!n) { 4784 PetscDrawAxis axis; 4785 4786 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4787 ierr = PetscDrawAxisSetLabels(axis,"Nonlinear iterations as function of time","Time","SNES Iterations");CHKERRQ(ierr); 4788 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4789 4790 ctx->snes_its = 0; 4791 } 4792 ierr = TSGetSNESIterations(ts,&its);CHKERRQ(ierr); 4793 y = its - ctx->snes_its; 4794 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4795 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 4796 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4797 } 4798 ctx->snes_its = its; 4799 PetscFunctionReturn(0); 4800 } 4801 4802 #undef __FUNCT__ 4803 #define __FUNCT__ "TSMonitorLGKSPIterations" 4804 PetscErrorCode TSMonitorLGKSPIterations(TS ts,PetscInt n,PetscReal ptime,Vec v,void *monctx) 4805 { 4806 TSMonitorLGCtx ctx = (TSMonitorLGCtx) monctx; 4807 PetscReal x = ptime,y; 4808 PetscErrorCode ierr; 4809 PetscInt its; 4810 4811 PetscFunctionBegin; 4812 if (!n) { 4813 PetscDrawAxis axis; 4814 4815 ierr = PetscDrawLGGetAxis(ctx->lg,&axis);CHKERRQ(ierr); 4816 ierr = PetscDrawAxisSetLabels(axis,"Linear iterations as function of time","Time","KSP Iterations");CHKERRQ(ierr); 4817 ierr = PetscDrawLGReset(ctx->lg);CHKERRQ(ierr); 4818 4819 ctx->ksp_its = 0; 4820 } 4821 ierr = TSGetKSPIterations(ts,&its);CHKERRQ(ierr); 4822 y = its - ctx->ksp_its; 4823 ierr = PetscDrawLGAddPoint(ctx->lg,&x,&y);CHKERRQ(ierr); 4824 if (((ctx->howoften > 0) && (!(n % ctx->howoften)) && (n > -1)) || ((ctx->howoften == -1) && (n == -1))) { 4825 ierr = PetscDrawLGDraw(ctx->lg);CHKERRQ(ierr); 4826 } 4827 ctx->ksp_its = its; 4828 PetscFunctionReturn(0); 4829 } 4830 4831 #undef __FUNCT__ 4832 #define __FUNCT__ "TSComputeLinearStability" 4833 /*@ 4834 TSComputeLinearStability - computes the linear stability function at a point 4835 4836 Collective on TS and Vec 4837 4838 Input Parameters: 4839 + ts - the TS context 4840 - xr,xi - real and imaginary part of input arguments 4841 4842 Output Parameters: 4843 . yr,yi - real and imaginary part of function value 4844 4845 Level: developer 4846 4847 .keywords: TS, compute 4848 4849 .seealso: TSSetRHSFunction(), TSComputeIFunction() 4850 @*/ 4851 PetscErrorCode TSComputeLinearStability(TS ts,PetscReal xr,PetscReal xi,PetscReal *yr,PetscReal *yi) 4852 { 4853 PetscErrorCode ierr; 4854 4855 PetscFunctionBegin; 4856 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 4857 if (!ts->ops->linearstability) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"Linearized stability function not provided for this method"); 4858 ierr = (*ts->ops->linearstability)(ts,xr,xi,yr,yi);CHKERRQ(ierr); 4859 PetscFunctionReturn(0); 4860 } 4861 4862 #undef __FUNCT__ 4863 #define __FUNCT__ "TSRollBack" 4864 /*@ 4865 TSRollBack - Rolls back one time step 4866 4867 Collective on TS 4868 4869 Input Parameter: 4870 . ts - the TS context obtained from TSCreate() 4871 4872 Level: advanced 4873 4874 .keywords: TS, timestep, rollback 4875 4876 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate() 4877 @*/ 4878 PetscErrorCode TSRollBack(TS ts) 4879 { 4880 PetscErrorCode ierr; 4881 4882 PetscFunctionBegin; 4883 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 4884 4885 if (!ts->ops->rollback) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSRollBack not implemented for type '%s'",((PetscObject)ts)->type_name); 4886 ierr = (*ts->ops->rollback)(ts);CHKERRQ(ierr); 4887 ts->time_step = ts->ptime - ts->ptime_prev; 4888 ts->ptime = ts->ptime_prev; 4889 PetscFunctionReturn(0); 4890 } 4891 4892 #undef __FUNCT__ 4893 #define __FUNCT__ "TSSetEventMonitor" 4894 /*@C 4895 TSSetEventMonitor - Sets a monitoring function used for detecting events 4896 4897 Logically Collective on TS 4898 4899 Input Parameters: 4900 + ts - the TS context obtained from TSCreate() 4901 . nevents - number of events 4902 . eventmonitor - event monitoring routine 4903 . postevent - [optional] post-event function 4904 - mectx - [optional] user-defined context for private data for the 4905 event monitor and post event routine (use NULL if no 4906 context is desired) 4907 4908 Calling sequence of eventmonitor: 4909 PetscErrorCode EventMonitor(TS ts,PetscReal t,Vec U,PetscScalar *fvalue,PetscInt *direction,PetscBool *terminate,void* mectx) 4910 4911 Input Parameters: 4912 + ts - the TS context 4913 . t - current time 4914 . U - current iterate 4915 - ctx - [optional] context passed with eventmonitor 4916 4917 Output parameters: 4918 + fvalue - function value of events at time t 4919 . direction - direction of zero crossing to be detected. -1 => Zero crossing in negative direction, 4920 +1 => Zero crossing in positive direction, 0 => both ways 4921 - terminate - terminate time stepping after event is detected. 4922 4923 Calling sequence of postevent: 4924 PetscErrorCode PostEvent(TS ts,PetscInt nevents_zero, PetscInt events_zero, PetscReal t,Vec U,void* ctx) 4925 4926 Input Parameters: 4927 + ts - the TS context 4928 . nevents_zero - number of events whose event function is zero 4929 . events_zero - indices of events which have reached zero 4930 . t - current time 4931 . U - current solution 4932 - ctx - the context passed with eventmonitor 4933 4934 Level: intermediate 4935 4936 .keywords: TS, event, set, monitor 4937 4938 .seealso: TSCreate(), TSSetTimeStep(), TSSetConvergedReason() 4939 @*/ 4940 PetscErrorCode TSSetEventMonitor(TS ts,PetscInt nevents,PetscErrorCode (*eventmonitor)(TS,PetscReal,Vec,PetscScalar*,PetscInt*,PetscBool*,void*),PetscErrorCode (*postevent)(TS,PetscInt,PetscInt[],PetscReal,Vec,void*),void *mectx) 4941 { 4942 PetscErrorCode ierr; 4943 PetscReal t; 4944 Vec U; 4945 TSEvent event; 4946 4947 PetscFunctionBegin; 4948 ierr = PetscNew(struct _p_TSEvent,&ts->event);CHKERRQ(ierr); 4949 event = ts->event; 4950 ierr = PetscMalloc(nevents*sizeof(PetscScalar),&event->fvalue);CHKERRQ(ierr); 4951 ierr = PetscMalloc(nevents*sizeof(PetscScalar),&event->fvalue_prev);CHKERRQ(ierr); 4952 ierr = PetscMalloc(nevents*sizeof(PetscBool),&event->terminate);CHKERRQ(ierr); 4953 ierr = PetscMalloc(nevents*sizeof(PetscInt),&event->direction);CHKERRQ(ierr); 4954 ierr = PetscMalloc(nevents*sizeof(PetscInt),&event->events_zero);CHKERRQ(ierr); 4955 event->monitor = eventmonitor; 4956 event->postevent = postevent; 4957 event->monitorcontext = (void*)mectx; 4958 event->nevents = nevents; 4959 4960 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 4961 ierr = TSGetTimeStep(ts,&event->initial_timestep);CHKERRQ(ierr); 4962 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 4963 event->ptime_prev = t; 4964 ierr = (*event->monitor)(ts,t,U,event->fvalue_prev,event->direction,event->terminate,mectx);CHKERRQ(ierr); 4965 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"TS Event options","");CHKERRQ(ierr); 4966 { 4967 event->tol = 1.0e-6; 4968 ierr = PetscOptionsReal("-ts_event_tol","","",event->tol,&event->tol,NULL);CHKERRQ(ierr); 4969 } 4970 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4971 PetscFunctionReturn(0); 4972 } 4973 4974 #undef __FUNCT__ 4975 #define __FUNCT__ "TSPostEvent" 4976 PetscErrorCode TSPostEvent(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,void *ctx) 4977 { 4978 PetscErrorCode ierr; 4979 TSEvent event=ts->event; 4980 PetscBool terminate=PETSC_FALSE; 4981 PetscInt i; 4982 PetscBool ts_terminate; 4983 4984 PetscFunctionBegin; 4985 if (event->postevent) { 4986 ierr = (*event->postevent)(ts,nevents_zero,events_zero,t,U,ctx);CHKERRQ(ierr); 4987 } 4988 for(i = 0; i < nevents_zero;i++) { 4989 terminate = terminate || event->terminate[events_zero[i]]; 4990 } 4991 ierr = MPI_Allreduce(&terminate,&ts_terminate,1,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 4992 if (terminate) { 4993 ierr = TSSetConvergedReason(ts,TS_CONVERGED_EVENT);CHKERRQ(ierr); 4994 event->status = TSEVENT_NONE; 4995 } else { 4996 event->status = TSEVENT_RESET_NEXTSTEP; 4997 } 4998 PetscFunctionReturn(0); 4999 } 5000 5001 #undef __FUNCT__ 5002 #define __FUNCT__ "TSEventMonitorDestroy" 5003 PetscErrorCode TSEventMonitorDestroy(TSEvent *event) 5004 { 5005 PetscErrorCode ierr; 5006 5007 PetscFunctionBegin; 5008 ierr = PetscFree((*event)->fvalue);CHKERRQ(ierr); 5009 ierr = PetscFree((*event)->fvalue_prev);CHKERRQ(ierr); 5010 ierr = PetscFree((*event)->terminate);CHKERRQ(ierr); 5011 ierr = PetscFree((*event)->direction);CHKERRQ(ierr); 5012 ierr = PetscFree((*event)->events_zero);CHKERRQ(ierr); 5013 ierr = PetscFree(*event);CHKERRQ(ierr); 5014 *event = NULL; 5015 PetscFunctionReturn(0); 5016 } 5017 5018 #undef __FUNCT__ 5019 #define __FUNCT__ "TSEventMonitor" 5020 PetscErrorCode TSEventMonitor(TS ts) 5021 { 5022 PetscErrorCode ierr; 5023 TSEvent event=ts->event; 5024 PetscReal t; 5025 Vec U; 5026 PetscInt i; 5027 PetscReal dt; 5028 PetscInt status=event->status; 5029 PetscInt rollback=0,in[2],out[2]; 5030 5031 PetscFunctionBegin; 5032 5033 ierr = TSGetTime(ts,&t);CHKERRQ(ierr); 5034 ierr = TSGetSolution(ts,&U);CHKERRQ(ierr); 5035 5036 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 5037 if (event->status == TSEVENT_RESET_NEXTSTEP) { 5038 /* Take initial time step */ 5039 dt = event->initial_timestep; 5040 ts->time_step = dt; 5041 event->status = TSEVENT_NONE; 5042 } 5043 5044 if (event->status == TSEVENT_NONE) { 5045 event->tstepend = t; 5046 } 5047 5048 event->nevents_zero = 0; 5049 5050 ierr = (*event->monitor)(ts,t,U,event->fvalue,event->direction,event->terminate,event->monitorcontext);CHKERRQ(ierr); 5051 if (event->status != TSEVENT_NONE) { 5052 for (i=0; i < event->nevents; i++) { 5053 if (PetscAbs(event->fvalue[i]) < event->tol) { 5054 event->status = TSEVENT_ZERO; 5055 event->events_zero[event->nevents_zero++] = i; 5056 } 5057 } 5058 } 5059 5060 status = event->status; 5061 ierr = MPI_Allreduce(&status,&event->status,1,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 5062 5063 if (event->status == TSEVENT_ZERO) { 5064 dt = event->tstepend-t; 5065 ts->time_step = dt; 5066 ierr = TSPostEvent(ts,event->nevents_zero,event->events_zero,t,U,event->monitorcontext);CHKERRQ(ierr); 5067 for (i = 0; i < event->nevents; i++) { 5068 event->fvalue_prev[i] = event->fvalue[i]; 5069 } 5070 event->ptime_prev = t; 5071 PetscFunctionReturn(0); 5072 } 5073 5074 for (i = 0; i < event->nevents; i++) { 5075 if ((event->direction[i] < 0 && PetscSign(event->fvalue[i]) <= 0 && PetscSign(event->fvalue_prev[i]) >= 0) || \ 5076 (event->direction[i] > 0 && PetscSign(event->fvalue[i]) >= 0 && PetscSign(event->fvalue_prev[i]) <= 0) || \ 5077 (event->direction[i] == 0 && PetscSign(event->fvalue[i])*PetscSign(event->fvalue_prev[i]) <= 0)) { 5078 5079 event->status = TSEVENT_LOCATED_INTERVAL; 5080 rollback = 1; 5081 /* Compute linearly interpolated new time step */ 5082 dt = PetscMin(dt,-event->fvalue_prev[i]*(t - event->ptime_prev)/(event->fvalue[i] - event->fvalue_prev[i])); 5083 } 5084 } 5085 in[0] = event->status; 5086 in[1] = rollback; 5087 ierr = MPI_Allreduce(in,out,2,MPIU_INT,MPI_MAX,((PetscObject)ts)->comm);CHKERRQ(ierr); 5088 5089 rollback = out[1]; 5090 if (rollback) { 5091 event->status = TSEVENT_LOCATED_INTERVAL; 5092 } 5093 5094 if (event->status == TSEVENT_LOCATED_INTERVAL) { 5095 ierr = TSRollBack(ts);CHKERRQ(ierr); 5096 event->status = TSEVENT_PROCESSING; 5097 } else { 5098 for (i = 0; i < event->nevents; i++) { 5099 event->fvalue_prev[i] = event->fvalue[i]; 5100 } 5101 event->ptime_prev = t; 5102 if (event->status == TSEVENT_PROCESSING) { 5103 dt = event->tstepend - event->ptime_prev; 5104 } 5105 } 5106 PetscReal time_step; 5107 ierr = MPI_Allreduce(&dt,&time_step,1,MPIU_REAL,MPI_MIN,((PetscObject)ts)->comm);CHKERRQ(ierr); 5108 ts->time_step = time_step; 5109 PetscFunctionReturn(0); 5110 } 5111 5112