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 if (A != B) { 575 ierr = MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&ts->Brhs);CHKERRQ(ierr); 576 } else { 577 ts->Brhs = ts->Arhs; 578 ierr = PetscObjectReference((PetscObject)ts->Arhs);CHKERRQ(ierr); 579 } 580 } 581 *Brhs = ts->Brhs; 582 } 583 PetscFunctionReturn(0); 584 } 585 586 #undef __FUNCT__ 587 #define __FUNCT__ "TSComputeIFunction" 588 /*@ 589 TSComputeIFunction - Evaluates the DAE residual written in implicit form F(t,U,Udot)=0 590 591 Collective on TS and Vec 592 593 Input Parameters: 594 + ts - the TS context 595 . t - current time 596 . U - state vector 597 . Udot - time derivative of state vector 598 - imex - flag indicates if the method is IMEX so that the RHSFunction should be kept separate 599 600 Output Parameter: 601 . Y - right hand side 602 603 Note: 604 Most users should not need to explicitly call this routine, as it 605 is used internally within the nonlinear solvers. 606 607 If the user did did not write their equations in implicit form, this 608 function recasts them in implicit form. 609 610 Level: developer 611 612 .keywords: TS, compute 613 614 .seealso: TSSetIFunction(), TSComputeRHSFunction() 615 @*/ 616 PetscErrorCode TSComputeIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec Y,PetscBool imex) 617 { 618 PetscErrorCode ierr; 619 TSIFunction ifunction; 620 TSRHSFunction rhsfunction; 621 void *ctx; 622 DM dm; 623 624 PetscFunctionBegin; 625 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 626 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 627 PetscValidHeaderSpecific(Udot,VEC_CLASSID,4); 628 PetscValidHeaderSpecific(Y,VEC_CLASSID,5); 629 630 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 631 ierr = DMTSGetIFunction(dm,&ifunction,&ctx);CHKERRQ(ierr); 632 ierr = DMTSGetRHSFunction(dm,&rhsfunction,NULL);CHKERRQ(ierr); 633 634 if (!rhsfunction && !ifunction) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSFunction() and / or TSSetIFunction()"); 635 636 ierr = PetscLogEventBegin(TS_FunctionEval,ts,U,Udot,Y);CHKERRQ(ierr); 637 if (ifunction) { 638 PetscStackPush("TS user implicit function"); 639 ierr = (*ifunction)(ts,t,U,Udot,Y,ctx);CHKERRQ(ierr); 640 PetscStackPop; 641 } 642 if (imex) { 643 if (!ifunction) { 644 ierr = VecCopy(Udot,Y);CHKERRQ(ierr); 645 } 646 } else if (rhsfunction) { 647 if (ifunction) { 648 Vec Frhs; 649 ierr = TSGetRHSVec_Private(ts,&Frhs);CHKERRQ(ierr); 650 ierr = TSComputeRHSFunction(ts,t,U,Frhs);CHKERRQ(ierr); 651 ierr = VecAXPY(Y,-1,Frhs);CHKERRQ(ierr); 652 } else { 653 ierr = TSComputeRHSFunction(ts,t,U,Y);CHKERRQ(ierr); 654 ierr = VecAYPX(Y,-1,Udot);CHKERRQ(ierr); 655 } 656 } 657 ierr = PetscLogEventEnd(TS_FunctionEval,ts,U,Udot,Y);CHKERRQ(ierr); 658 PetscFunctionReturn(0); 659 } 660 661 #undef __FUNCT__ 662 #define __FUNCT__ "TSComputeIJacobian" 663 /*@ 664 TSComputeIJacobian - Evaluates the Jacobian of the DAE 665 666 Collective on TS and Vec 667 668 Input 669 Input Parameters: 670 + ts - the TS context 671 . t - current timestep 672 . U - state vector 673 . Udot - time derivative of state vector 674 . shift - shift to apply, see note below 675 - imex - flag indicates if the method is IMEX so that the RHSJacobian should be kept separate 676 677 Output Parameters: 678 + A - Jacobian matrix 679 . B - optional preconditioning matrix 680 - flag - flag indicating matrix structure 681 682 Notes: 683 If F(t,U,Udot)=0 is the DAE, the required Jacobian is 684 685 dF/dU + shift*dF/dUdot 686 687 Most users should not need to explicitly call this routine, as it 688 is used internally within the nonlinear solvers. 689 690 Level: developer 691 692 .keywords: TS, compute, Jacobian, matrix 693 694 .seealso: TSSetIJacobian() 695 @*/ 696 PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,PetscBool imex) 697 { 698 PetscErrorCode ierr; 699 TSIJacobian ijacobian; 700 TSRHSJacobian rhsjacobian; 701 DM dm; 702 void *ctx; 703 704 PetscFunctionBegin; 705 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 706 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 707 PetscValidHeaderSpecific(Udot,VEC_CLASSID,4); 708 PetscValidPointer(A,6); 709 PetscValidHeaderSpecific(*A,MAT_CLASSID,6); 710 PetscValidPointer(B,7); 711 PetscValidHeaderSpecific(*B,MAT_CLASSID,7); 712 PetscValidPointer(flg,8); 713 714 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 715 ierr = DMTSGetIJacobian(dm,&ijacobian,&ctx);CHKERRQ(ierr); 716 ierr = DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);CHKERRQ(ierr); 717 718 if (!rhsjacobian && !ijacobian) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()"); 719 720 *flg = SAME_NONZERO_PATTERN; /* In case we're solving a linear problem in which case it wouldn't get initialized below. */ 721 ierr = PetscLogEventBegin(TS_JacobianEval,ts,U,*A,*B);CHKERRQ(ierr); 722 if (ijacobian) { 723 *flg = DIFFERENT_NONZERO_PATTERN; 724 PetscStackPush("TS user implicit Jacobian"); 725 ierr = (*ijacobian)(ts,t,U,Udot,shift,A,B,flg,ctx);CHKERRQ(ierr); 726 PetscStackPop; 727 /* make sure user returned a correct Jacobian and preconditioner */ 728 PetscValidHeaderSpecific(*A,MAT_CLASSID,4); 729 PetscValidHeaderSpecific(*B,MAT_CLASSID,5); 730 } 731 if (imex) { 732 if (!ijacobian) { /* system was written as Udot = G(t,U) */ 733 ierr = MatZeroEntries(*A);CHKERRQ(ierr); 734 ierr = MatShift(*A,shift);CHKERRQ(ierr); 735 if (*A != *B) { 736 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 737 ierr = MatShift(*B,shift);CHKERRQ(ierr); 738 } 739 *flg = SAME_PRECONDITIONER; 740 } 741 } else { 742 Mat Arhs = NULL,Brhs = NULL; 743 MatStructure flg2; 744 if (rhsjacobian) { 745 ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr); 746 ierr = TSComputeRHSJacobian(ts,t,U,&Arhs,&Brhs,&flg2);CHKERRQ(ierr); 747 } 748 if (Arhs == *A) { /* No IJacobian, so we only have the RHS matrix */ 749 ts->rhsjacobian.scale = -1; 750 ts->rhsjacobian.shift = shift; 751 ierr = MatScale(*A,-1);CHKERRQ(ierr); 752 ierr = MatShift(*A,shift);CHKERRQ(ierr); 753 if (*A != *B) { 754 ierr = MatScale(*B,-1);CHKERRQ(ierr); 755 ierr = MatShift(*B,shift);CHKERRQ(ierr); 756 } 757 } else if (Arhs) { /* Both IJacobian and RHSJacobian */ 758 MatStructure axpy = DIFFERENT_NONZERO_PATTERN; 759 if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */ 760 ierr = MatZeroEntries(*A);CHKERRQ(ierr); 761 ierr = MatShift(*A,shift);CHKERRQ(ierr); 762 if (*A != *B) { 763 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 764 ierr = MatShift(*B,shift);CHKERRQ(ierr); 765 } 766 } 767 ierr = MatAXPY(*A,-1,Arhs,axpy);CHKERRQ(ierr); 768 if (*A != *B) { 769 ierr = MatAXPY(*B,-1,Brhs,axpy);CHKERRQ(ierr); 770 } 771 *flg = PetscMin(*flg,flg2); 772 } 773 } 774 775 ierr = PetscLogEventEnd(TS_JacobianEval,ts,U,*A,*B);CHKERRQ(ierr); 776 PetscFunctionReturn(0); 777 } 778 779 #undef __FUNCT__ 780 #define __FUNCT__ "TSSetRHSFunction" 781 /*@C 782 TSSetRHSFunction - Sets the routine for evaluating the function, 783 where U_t = G(t,u). 784 785 Logically Collective on TS 786 787 Input Parameters: 788 + ts - the TS context obtained from TSCreate() 789 . r - vector to put the computed right hand side (or NULL to have it created) 790 . f - routine for evaluating the right-hand-side function 791 - ctx - [optional] user-defined context for private data for the 792 function evaluation routine (may be NULL) 793 794 Calling sequence of func: 795 $ func (TS ts,PetscReal t,Vec u,Vec F,void *ctx); 796 797 + t - current timestep 798 . u - input vector 799 . F - function vector 800 - ctx - [optional] user-defined function context 801 802 Level: beginner 803 804 .keywords: TS, timestep, set, right-hand-side, function 805 806 .seealso: TSSetRHSJacobian(), TSSetIJacobian() 807 @*/ 808 PetscErrorCode TSSetRHSFunction(TS ts,Vec r,PetscErrorCode (*f)(TS,PetscReal,Vec,Vec,void*),void *ctx) 809 { 810 PetscErrorCode ierr; 811 SNES snes; 812 Vec ralloc = NULL; 813 DM dm; 814 815 PetscFunctionBegin; 816 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 817 if (r) PetscValidHeaderSpecific(r,VEC_CLASSID,2); 818 819 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 820 ierr = DMTSSetRHSFunction(dm,f,ctx);CHKERRQ(ierr); 821 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 822 if (!r && !ts->dm && ts->vec_sol) { 823 ierr = VecDuplicate(ts->vec_sol,&ralloc);CHKERRQ(ierr); 824 r = ralloc; 825 } 826 ierr = SNESSetFunction(snes,r,SNESTSFormFunction,ts);CHKERRQ(ierr); 827 ierr = VecDestroy(&ralloc);CHKERRQ(ierr); 828 PetscFunctionReturn(0); 829 } 830 831 #undef __FUNCT__ 832 #define __FUNCT__ "TSSetSolutionFunction" 833 /*@C 834 TSSetSolutionFunction - Provide a function that computes the solution of the ODE or DAE 835 836 Logically Collective on TS 837 838 Input Parameters: 839 + ts - the TS context obtained from TSCreate() 840 . f - routine for evaluating the solution 841 - ctx - [optional] user-defined context for private data for the 842 function evaluation routine (may be NULL) 843 844 Calling sequence of func: 845 $ func (TS ts,PetscReal t,Vec u,void *ctx); 846 847 + t - current timestep 848 . u - output vector 849 - ctx - [optional] user-defined function context 850 851 Notes: 852 This routine is used for testing accuracy of time integration schemes when you already know the solution. 853 If analytic solutions are not known for your system, consider using the Method of Manufactured Solutions to 854 create closed-form solutions with non-physical forcing terms. 855 856 For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history. 857 858 Level: beginner 859 860 .keywords: TS, timestep, set, right-hand-side, function 861 862 .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetForcingFunction() 863 @*/ 864 PetscErrorCode TSSetSolutionFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx) 865 { 866 PetscErrorCode ierr; 867 DM dm; 868 869 PetscFunctionBegin; 870 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 871 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 872 ierr = DMTSSetSolutionFunction(dm,f,ctx);CHKERRQ(ierr); 873 PetscFunctionReturn(0); 874 } 875 876 #undef __FUNCT__ 877 #define __FUNCT__ "TSSetForcingFunction" 878 /*@C 879 TSSetForcingFunction - Provide a function that computes a forcing term for a ODE or PDE 880 881 Logically Collective on TS 882 883 Input Parameters: 884 + ts - the TS context obtained from TSCreate() 885 . f - routine for evaluating the forcing function 886 - ctx - [optional] user-defined context for private data for the 887 function evaluation routine (may be NULL) 888 889 Calling sequence of func: 890 $ func (TS ts,PetscReal t,Vec u,void *ctx); 891 892 + t - current timestep 893 . u - output vector 894 - ctx - [optional] user-defined function context 895 896 Notes: 897 This routine is useful for testing accuracy of time integration schemes when using the Method of Manufactured Solutions to 898 create closed-form solutions with a non-physical forcing term. 899 900 For low-dimensional problems solved in serial, such as small discrete systems, TSMonitorLGError() can be used to monitor the error history. 901 902 Level: beginner 903 904 .keywords: TS, timestep, set, right-hand-side, function 905 906 .seealso: TSSetRHSJacobian(), TSSetIJacobian(), TSComputeSolutionFunction(), TSSetSolutionFunction() 907 @*/ 908 PetscErrorCode TSSetForcingFunction(TS ts,PetscErrorCode (*f)(TS,PetscReal,Vec,void*),void *ctx) 909 { 910 PetscErrorCode ierr; 911 DM dm; 912 913 PetscFunctionBegin; 914 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 915 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 916 ierr = DMTSSetForcingFunction(dm,f,ctx);CHKERRQ(ierr); 917 PetscFunctionReturn(0); 918 } 919 920 #undef __FUNCT__ 921 #define __FUNCT__ "TSSetRHSJacobian" 922 /*@C 923 TSSetRHSJacobian - Sets the function to compute the Jacobian of F, 924 where U_t = G(U,t), as well as the location to store the matrix. 925 926 Logically Collective on TS 927 928 Input Parameters: 929 + ts - the TS context obtained from TSCreate() 930 . Amat - (approximate) Jacobian matrix 931 . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat) 932 . f - the Jacobian evaluation routine 933 - ctx - [optional] user-defined context for private data for the 934 Jacobian evaluation routine (may be NULL) 935 936 Calling sequence of func: 937 $ func (TS ts,PetscReal t,Vec u,Mat *A,Mat *B,MatStructure *flag,void *ctx); 938 939 + t - current timestep 940 . u - input vector 941 . Amat - (approximate) Jacobian matrix 942 . Pmat - matrix from which preconditioner is to be constructed (usually the same as Amat) 943 . flag - flag indicating information about the preconditioner matrix 944 structure (same as flag in KSPSetOperators()) 945 - ctx - [optional] user-defined context for matrix evaluation routine 946 947 Notes: 948 See KSPSetOperators() for important information about setting the flag 949 output parameter in the routine func(). Be sure to read this information! 950 951 The routine func() takes Mat * as the matrix arguments rather than Mat. 952 This allows the matrix evaluation routine to replace A and/or B with a 953 completely new matrix structure (not just different matrix elements) 954 when appropriate, for instance, if the nonzero structure is changing 955 throughout the global iterations. 956 957 Level: beginner 958 959 .keywords: TS, timestep, set, right-hand-side, Jacobian 960 961 .seealso: SNESComputeJacobianDefaultColor(), TSSetRHSFunction(), TSRHSJacobianSetReuse() 962 963 @*/ 964 PetscErrorCode TSSetRHSJacobian(TS ts,Mat Amat,Mat Pmat,TSRHSJacobian f,void *ctx) 965 { 966 PetscErrorCode ierr; 967 SNES snes; 968 DM dm; 969 TSIJacobian ijacobian; 970 971 PetscFunctionBegin; 972 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 973 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 974 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 975 if (Amat) PetscCheckSameComm(ts,1,Amat,2); 976 if (Pmat) PetscCheckSameComm(ts,1,Pmat,3); 977 978 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 979 ierr = DMTSSetRHSJacobian(dm,f,ctx);CHKERRQ(ierr); 980 if (f == TSComputeRHSJacobianConstant) { 981 /* Handle this case automatically for the user; otherwise user should call themselves. */ 982 ierr = TSRHSJacobianSetReuse(ts,PETSC_TRUE);CHKERRQ(ierr); 983 } 984 ierr = DMTSGetIJacobian(dm,&ijacobian,NULL);CHKERRQ(ierr); 985 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 986 if (!ijacobian) { 987 ierr = SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);CHKERRQ(ierr); 988 } 989 if (Amat) { 990 ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr); 991 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 992 993 ts->Arhs = Amat; 994 } 995 if (Pmat) { 996 ierr = PetscObjectReference((PetscObject)Pmat);CHKERRQ(ierr); 997 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 998 999 ts->Brhs = Pmat; 1000 } 1001 PetscFunctionReturn(0); 1002 } 1003 1004 1005 #undef __FUNCT__ 1006 #define __FUNCT__ "TSSetIFunction" 1007 /*@C 1008 TSSetIFunction - Set the function to compute F(t,U,U_t) where F() = 0 is the DAE to be solved. 1009 1010 Logically Collective on TS 1011 1012 Input Parameters: 1013 + ts - the TS context obtained from TSCreate() 1014 . r - vector to hold the residual (or NULL to have it created internally) 1015 . f - the function evaluation routine 1016 - ctx - user-defined context for private data for the function evaluation routine (may be NULL) 1017 1018 Calling sequence of f: 1019 $ f(TS ts,PetscReal t,Vec u,Vec u_t,Vec F,ctx); 1020 1021 + t - time at step/stage being solved 1022 . u - state vector 1023 . u_t - time derivative of state vector 1024 . F - function vector 1025 - ctx - [optional] user-defined context for matrix evaluation routine 1026 1027 Important: 1028 The user MUST call either this routine, TSSetRHSFunction(). This routine must be used when not solving an ODE, for example a DAE. 1029 1030 Level: beginner 1031 1032 .keywords: TS, timestep, set, DAE, Jacobian 1033 1034 .seealso: TSSetRHSJacobian(), TSSetRHSFunction(), TSSetIJacobian() 1035 @*/ 1036 PetscErrorCode TSSetIFunction(TS ts,Vec res,TSIFunction f,void *ctx) 1037 { 1038 PetscErrorCode ierr; 1039 SNES snes; 1040 Vec resalloc = NULL; 1041 DM dm; 1042 1043 PetscFunctionBegin; 1044 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1045 if (res) PetscValidHeaderSpecific(res,VEC_CLASSID,2); 1046 1047 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1048 ierr = DMTSSetIFunction(dm,f,ctx);CHKERRQ(ierr); 1049 1050 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1051 if (!res && !ts->dm && ts->vec_sol) { 1052 ierr = VecDuplicate(ts->vec_sol,&resalloc);CHKERRQ(ierr); 1053 res = resalloc; 1054 } 1055 ierr = SNESSetFunction(snes,res,SNESTSFormFunction,ts);CHKERRQ(ierr); 1056 ierr = VecDestroy(&resalloc);CHKERRQ(ierr); 1057 PetscFunctionReturn(0); 1058 } 1059 1060 #undef __FUNCT__ 1061 #define __FUNCT__ "TSGetIFunction" 1062 /*@C 1063 TSGetIFunction - Returns the vector where the implicit residual is stored and the function/contex to compute it. 1064 1065 Not Collective 1066 1067 Input Parameter: 1068 . ts - the TS context 1069 1070 Output Parameter: 1071 + r - vector to hold residual (or NULL) 1072 . func - the function to compute residual (or NULL) 1073 - ctx - the function context (or NULL) 1074 1075 Level: advanced 1076 1077 .keywords: TS, nonlinear, get, function 1078 1079 .seealso: TSSetIFunction(), SNESGetFunction() 1080 @*/ 1081 PetscErrorCode TSGetIFunction(TS ts,Vec *r,TSIFunction *func,void **ctx) 1082 { 1083 PetscErrorCode ierr; 1084 SNES snes; 1085 DM dm; 1086 1087 PetscFunctionBegin; 1088 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1089 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1090 ierr = SNESGetFunction(snes,r,NULL,NULL);CHKERRQ(ierr); 1091 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1092 ierr = DMTSGetIFunction(dm,func,ctx);CHKERRQ(ierr); 1093 PetscFunctionReturn(0); 1094 } 1095 1096 #undef __FUNCT__ 1097 #define __FUNCT__ "TSGetRHSFunction" 1098 /*@C 1099 TSGetRHSFunction - Returns the vector where the right hand side is stored and the function/context to compute it. 1100 1101 Not Collective 1102 1103 Input Parameter: 1104 . ts - the TS context 1105 1106 Output Parameter: 1107 + r - vector to hold computed right hand side (or NULL) 1108 . func - the function to compute right hand side (or NULL) 1109 - ctx - the function context (or NULL) 1110 1111 Level: advanced 1112 1113 .keywords: TS, nonlinear, get, function 1114 1115 .seealso: TSSetRhsfunction(), SNESGetFunction() 1116 @*/ 1117 PetscErrorCode TSGetRHSFunction(TS ts,Vec *r,TSRHSFunction *func,void **ctx) 1118 { 1119 PetscErrorCode ierr; 1120 SNES snes; 1121 DM dm; 1122 1123 PetscFunctionBegin; 1124 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1125 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1126 ierr = SNESGetFunction(snes,r,NULL,NULL);CHKERRQ(ierr); 1127 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1128 ierr = DMTSGetRHSFunction(dm,func,ctx);CHKERRQ(ierr); 1129 PetscFunctionReturn(0); 1130 } 1131 1132 #undef __FUNCT__ 1133 #define __FUNCT__ "TSSetIJacobian" 1134 /*@C 1135 TSSetIJacobian - Set the function to compute the matrix dF/dU + a*dF/dU_t where F(t,U,U_t) is the function 1136 you provided with TSSetIFunction(). 1137 1138 Logically Collective on TS 1139 1140 Input Parameters: 1141 + ts - the TS context obtained from TSCreate() 1142 . Amat - (approximate) Jacobian matrix 1143 . Pmat - matrix used to compute preconditioner (usually the same as Amat) 1144 . f - the Jacobian evaluation routine 1145 - ctx - user-defined context for private data for the Jacobian evaluation routine (may be NULL) 1146 1147 Calling sequence of f: 1148 $ f(TS ts,PetscReal t,Vec U,Vec U_t,PetscReal a,Mat *Amat,Mat *Pmat,MatStructure *flag,void *ctx); 1149 1150 + t - time at step/stage being solved 1151 . U - state vector 1152 . U_t - time derivative of state vector 1153 . a - shift 1154 . Amat - (approximate) Jacobian of F(t,U,W+a*U), equivalent to dF/dU + a*dF/dU_t 1155 . Pmat - matrix used for constructing preconditioner, usually the same as Amat 1156 . flag - flag indicating information about the preconditioner matrix 1157 structure (same as flag in KSPSetOperators()) 1158 - ctx - [optional] user-defined context for matrix evaluation routine 1159 1160 Notes: 1161 The matrices Amat and Pmat are exactly the matrices that are used by SNES for the nonlinear solve. 1162 1163 The matrix dF/dU + a*dF/dU_t you provide turns out to be 1164 the Jacobian of F(t,U,W+a*U) where F(t,U,U_t) = 0 is the DAE to be solved. 1165 The time integrator internally approximates U_t by W+a*U where the positive "shift" 1166 a and vector W depend on the integration method, step size, and past states. For example with 1167 the backward Euler method a = 1/dt and W = -a*U(previous timestep) so 1168 W + a*U = a*(U - U(previous timestep)) = (U - U(previous timestep))/dt 1169 1170 Level: beginner 1171 1172 .keywords: TS, timestep, DAE, Jacobian 1173 1174 .seealso: TSSetIFunction(), TSSetRHSJacobian(), SNESComputeJacobianDefaultColor(), SNESComputeJacobianDefault() 1175 1176 @*/ 1177 PetscErrorCode TSSetIJacobian(TS ts,Mat Amat,Mat Pmat,TSIJacobian f,void *ctx) 1178 { 1179 PetscErrorCode ierr; 1180 SNES snes; 1181 DM dm; 1182 1183 PetscFunctionBegin; 1184 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1185 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 1186 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 1187 if (Amat) PetscCheckSameComm(ts,1,Amat,2); 1188 if (Pmat) PetscCheckSameComm(ts,1,Pmat,3); 1189 1190 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1191 ierr = DMTSSetIJacobian(dm,f,ctx);CHKERRQ(ierr); 1192 1193 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1194 ierr = SNESSetJacobian(snes,Amat,Pmat,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1195 PetscFunctionReturn(0); 1196 } 1197 1198 #undef __FUNCT__ 1199 #define __FUNCT__ "TSRHSJacobianSetReuse" 1200 /*@ 1201 TSRHSJacobianSetReuse - restore RHS Jacobian before re-evaluating. Without this flag, TS will change the sign and 1202 shift the RHS Jacobian for a finite-time-step implicit solve, in which case the user function will need to recompute 1203 the entire Jacobian. The reuse flag must be set if the evaluation function will assume that the matrix entries have 1204 not been changed by the TS. 1205 1206 Logically Collective 1207 1208 Input Arguments: 1209 + ts - TS context obtained from TSCreate() 1210 - reuse - PETSC_TRUE if the RHS Jacobian 1211 1212 Level: intermediate 1213 1214 .seealso: TSSetRHSJacobian(), TSComputeRHSJacobianConstant() 1215 @*/ 1216 PetscErrorCode TSRHSJacobianSetReuse(TS ts,PetscBool reuse) 1217 { 1218 PetscFunctionBegin; 1219 ts->rhsjacobian.reuse = reuse; 1220 PetscFunctionReturn(0); 1221 } 1222 1223 #undef __FUNCT__ 1224 #define __FUNCT__ "TSLoad" 1225 /*@C 1226 TSLoad - Loads a KSP that has been stored in binary with KSPView(). 1227 1228 Collective on PetscViewer 1229 1230 Input Parameters: 1231 + newdm - the newly loaded TS, this needs to have been created with TSCreate() or 1232 some related function before a call to TSLoad(). 1233 - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() 1234 1235 Level: intermediate 1236 1237 Notes: 1238 The type is determined by the data in the file, any type set into the TS before this call is ignored. 1239 1240 Notes for advanced users: 1241 Most users should not need to know the details of the binary storage 1242 format, since TSLoad() and TSView() completely hide these details. 1243 But for anyone who's interested, the standard binary matrix storage 1244 format is 1245 .vb 1246 has not yet been determined 1247 .ve 1248 1249 .seealso: PetscViewerBinaryOpen(), TSView(), MatLoad(), VecLoad() 1250 @*/ 1251 PetscErrorCode TSLoad(TS ts, PetscViewer viewer) 1252 { 1253 PetscErrorCode ierr; 1254 PetscBool isbinary; 1255 PetscInt classid; 1256 char type[256]; 1257 DMTS sdm; 1258 DM dm; 1259 1260 PetscFunctionBegin; 1261 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1262 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1263 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1264 if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()"); 1265 1266 ierr = PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr); 1267 if (classid != TS_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_WRONG,"Not TS next in file"); 1268 ierr = PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);CHKERRQ(ierr); 1269 ierr = TSSetType(ts, type);CHKERRQ(ierr); 1270 if (ts->ops->load) { 1271 ierr = (*ts->ops->load)(ts,viewer);CHKERRQ(ierr); 1272 } 1273 ierr = DMCreate(PetscObjectComm((PetscObject)ts),&dm);CHKERRQ(ierr); 1274 ierr = DMLoad(dm,viewer);CHKERRQ(ierr); 1275 ierr = TSSetDM(ts,dm);CHKERRQ(ierr); 1276 ierr = DMCreateGlobalVector(ts->dm,&ts->vec_sol);CHKERRQ(ierr); 1277 ierr = VecLoad(ts->vec_sol,viewer);CHKERRQ(ierr); 1278 ierr = DMGetDMTS(ts->dm,&sdm);CHKERRQ(ierr); 1279 ierr = DMTSLoad(sdm,viewer);CHKERRQ(ierr); 1280 PetscFunctionReturn(0); 1281 } 1282 1283 #include <petscdraw.h> 1284 #if defined(PETSC_HAVE_AMS) 1285 #include <petscviewerams.h> 1286 #endif 1287 #undef __FUNCT__ 1288 #define __FUNCT__ "TSView" 1289 /*@C 1290 TSView - Prints the TS data structure. 1291 1292 Collective on TS 1293 1294 Input Parameters: 1295 + ts - the TS context obtained from TSCreate() 1296 - viewer - visualization context 1297 1298 Options Database Key: 1299 . -ts_view - calls TSView() at end of TSStep() 1300 1301 Notes: 1302 The available visualization contexts include 1303 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 1304 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 1305 output where only the first processor opens 1306 the file. All other processors send their 1307 data to the first processor to print. 1308 1309 The user can open an alternative visualization context with 1310 PetscViewerASCIIOpen() - output to a specified file. 1311 1312 Level: beginner 1313 1314 .keywords: TS, timestep, view 1315 1316 .seealso: PetscViewerASCIIOpen() 1317 @*/ 1318 PetscErrorCode TSView(TS ts,PetscViewer viewer) 1319 { 1320 PetscErrorCode ierr; 1321 TSType type; 1322 PetscBool iascii,isstring,isundials,isbinary,isdraw; 1323 DMTS sdm; 1324 #if defined(PETSC_HAVE_AMS) 1325 PetscBool isams; 1326 #endif 1327 1328 PetscFunctionBegin; 1329 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1330 if (!viewer) { 1331 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ts),&viewer);CHKERRQ(ierr); 1332 } 1333 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1334 PetscCheckSameComm(ts,1,viewer,2); 1335 1336 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1337 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 1338 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1339 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1340 #if defined(PETSC_HAVE_AMS) 1341 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERAMS,&isams);CHKERRQ(ierr); 1342 #endif 1343 if (iascii) { 1344 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ts,viewer,"TS Object");CHKERRQ(ierr); 1345 ierr = PetscViewerASCIIPrintf(viewer," maximum steps=%D\n",ts->max_steps);CHKERRQ(ierr); 1346 ierr = PetscViewerASCIIPrintf(viewer," maximum time=%G\n",ts->max_time);CHKERRQ(ierr); 1347 if (ts->problem_type == TS_NONLINEAR) { 1348 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solver iterations=%D\n",ts->snes_its);CHKERRQ(ierr); 1349 ierr = PetscViewerASCIIPrintf(viewer," total number of nonlinear solve failures=%D\n",ts->num_snes_failures);CHKERRQ(ierr); 1350 } 1351 ierr = PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",ts->ksp_its);CHKERRQ(ierr); 1352 ierr = PetscViewerASCIIPrintf(viewer," total number of rejected steps=%D\n",ts->reject);CHKERRQ(ierr); 1353 ierr = DMGetDMTS(ts->dm,&sdm);CHKERRQ(ierr); 1354 ierr = DMTSView(sdm,viewer);CHKERRQ(ierr); 1355 if (ts->ops->view) { 1356 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1357 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1358 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1359 } 1360 } else if (isstring) { 1361 ierr = TSGetType(ts,&type);CHKERRQ(ierr); 1362 ierr = PetscViewerStringSPrintf(viewer," %-7.7s",type);CHKERRQ(ierr); 1363 } else if (isbinary) { 1364 PetscInt classid = TS_FILE_CLASSID; 1365 MPI_Comm comm; 1366 PetscMPIInt rank; 1367 char type[256]; 1368 1369 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 1370 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1371 if (!rank) { 1372 ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1373 ierr = PetscStrncpy(type,((PetscObject)ts)->type_name,256);CHKERRQ(ierr); 1374 ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr); 1375 } 1376 if (ts->ops->view) { 1377 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1378 } 1379 ierr = DMView(ts->dm,viewer);CHKERRQ(ierr); 1380 ierr = VecView(ts->vec_sol,viewer);CHKERRQ(ierr); 1381 ierr = DMGetDMTS(ts->dm,&sdm);CHKERRQ(ierr); 1382 ierr = DMTSView(sdm,viewer);CHKERRQ(ierr); 1383 } else if (isdraw) { 1384 PetscDraw draw; 1385 char str[36]; 1386 PetscReal x,y,bottom,h; 1387 1388 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1389 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 1390 ierr = PetscStrcpy(str,"TS: ");CHKERRQ(ierr); 1391 ierr = PetscStrcat(str,((PetscObject)ts)->type_name);CHKERRQ(ierr); 1392 ierr = PetscDrawBoxedString(draw,x,y,PETSC_DRAW_BLACK,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 1393 bottom = y - h; 1394 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr); 1395 if (ts->ops->view) { 1396 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1397 } 1398 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 1399 #if defined(PETSC_HAVE_AMS) 1400 } else if (isams) { 1401 if (((PetscObject)ts)->amsmem == -1) { 1402 ierr = PetscObjectViewAMS((PetscObject)ts,viewer);CHKERRQ(ierr); 1403 PetscStackCallAMS(AMS_Memory_take_access,(((PetscObject)ts)->amsmem)); 1404 PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)ts)->amsmem,"time step",&ts->steps,1,AMS_INT,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF)); 1405 PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)ts)->amsmem,"time",&ts->ptime,1,AMS_DOUBLE,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF)); 1406 PetscStackCallAMS(AMS_Memory_grant_access,(((PetscObject)ts)->amsmem)); 1407 } 1408 if (ts->ops->view) { 1409 ierr = (*ts->ops->view)(ts,viewer);CHKERRQ(ierr); 1410 } 1411 #endif 1412 } 1413 1414 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1415 ierr = PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&isundials);CHKERRQ(ierr); 1416 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1417 PetscFunctionReturn(0); 1418 } 1419 1420 1421 #undef __FUNCT__ 1422 #define __FUNCT__ "TSSetApplicationContext" 1423 /*@ 1424 TSSetApplicationContext - Sets an optional user-defined context for 1425 the timesteppers. 1426 1427 Logically Collective on TS 1428 1429 Input Parameters: 1430 + ts - the TS context obtained from TSCreate() 1431 - usrP - optional user context 1432 1433 Level: intermediate 1434 1435 .keywords: TS, timestep, set, application, context 1436 1437 .seealso: TSGetApplicationContext() 1438 @*/ 1439 PetscErrorCode TSSetApplicationContext(TS ts,void *usrP) 1440 { 1441 PetscFunctionBegin; 1442 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1443 ts->user = usrP; 1444 PetscFunctionReturn(0); 1445 } 1446 1447 #undef __FUNCT__ 1448 #define __FUNCT__ "TSGetApplicationContext" 1449 /*@ 1450 TSGetApplicationContext - Gets the user-defined context for the 1451 timestepper. 1452 1453 Not Collective 1454 1455 Input Parameter: 1456 . ts - the TS context obtained from TSCreate() 1457 1458 Output Parameter: 1459 . usrP - user context 1460 1461 Level: intermediate 1462 1463 .keywords: TS, timestep, get, application, context 1464 1465 .seealso: TSSetApplicationContext() 1466 @*/ 1467 PetscErrorCode TSGetApplicationContext(TS ts,void *usrP) 1468 { 1469 PetscFunctionBegin; 1470 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1471 *(void**)usrP = ts->user; 1472 PetscFunctionReturn(0); 1473 } 1474 1475 #undef __FUNCT__ 1476 #define __FUNCT__ "TSGetTimeStepNumber" 1477 /*@ 1478 TSGetTimeStepNumber - Gets the number of time steps completed. 1479 1480 Not Collective 1481 1482 Input Parameter: 1483 . ts - the TS context obtained from TSCreate() 1484 1485 Output Parameter: 1486 . iter - number of steps completed so far 1487 1488 Level: intermediate 1489 1490 .keywords: TS, timestep, get, iteration, number 1491 .seealso: TSGetTime(), TSGetTimeStep(), TSSetPreStep(), TSSetPreStage(), TSSetPostStep() 1492 @*/ 1493 PetscErrorCode TSGetTimeStepNumber(TS ts,PetscInt *iter) 1494 { 1495 PetscFunctionBegin; 1496 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1497 PetscValidIntPointer(iter,2); 1498 *iter = ts->steps; 1499 PetscFunctionReturn(0); 1500 } 1501 1502 #undef __FUNCT__ 1503 #define __FUNCT__ "TSSetInitialTimeStep" 1504 /*@ 1505 TSSetInitialTimeStep - Sets the initial timestep to be used, 1506 as well as the initial time. 1507 1508 Logically Collective on TS 1509 1510 Input Parameters: 1511 + ts - the TS context obtained from TSCreate() 1512 . initial_time - the initial time 1513 - time_step - the size of the timestep 1514 1515 Level: intermediate 1516 1517 .seealso: TSSetTimeStep(), TSGetTimeStep() 1518 1519 .keywords: TS, set, initial, timestep 1520 @*/ 1521 PetscErrorCode TSSetInitialTimeStep(TS ts,PetscReal initial_time,PetscReal time_step) 1522 { 1523 PetscErrorCode ierr; 1524 1525 PetscFunctionBegin; 1526 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1527 ierr = TSSetTimeStep(ts,time_step);CHKERRQ(ierr); 1528 ierr = TSSetTime(ts,initial_time);CHKERRQ(ierr); 1529 PetscFunctionReturn(0); 1530 } 1531 1532 #undef __FUNCT__ 1533 #define __FUNCT__ "TSSetTimeStep" 1534 /*@ 1535 TSSetTimeStep - Allows one to reset the timestep at any time, 1536 useful for simple pseudo-timestepping codes. 1537 1538 Logically Collective on TS 1539 1540 Input Parameters: 1541 + ts - the TS context obtained from TSCreate() 1542 - time_step - the size of the timestep 1543 1544 Level: intermediate 1545 1546 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1547 1548 .keywords: TS, set, timestep 1549 @*/ 1550 PetscErrorCode TSSetTimeStep(TS ts,PetscReal time_step) 1551 { 1552 PetscFunctionBegin; 1553 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1554 PetscValidLogicalCollectiveReal(ts,time_step,2); 1555 ts->time_step = time_step; 1556 ts->time_step_orig = time_step; 1557 PetscFunctionReturn(0); 1558 } 1559 1560 #undef __FUNCT__ 1561 #define __FUNCT__ "TSSetExactFinalTime" 1562 /*@ 1563 TSSetExactFinalTime - Determines whether to adapt the final time step to 1564 match the exact final time, interpolate solution to the exact final time, 1565 or just return at the final time TS computed. 1566 1567 Logically Collective on TS 1568 1569 Input Parameter: 1570 + ts - the time-step context 1571 - eftopt - exact final time option 1572 1573 Level: beginner 1574 1575 .seealso: TSExactFinalTimeOption 1576 @*/ 1577 PetscErrorCode TSSetExactFinalTime(TS ts,TSExactFinalTimeOption eftopt) 1578 { 1579 PetscFunctionBegin; 1580 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1581 PetscValidLogicalCollectiveEnum(ts,eftopt,2); 1582 ts->exact_final_time = eftopt; 1583 PetscFunctionReturn(0); 1584 } 1585 1586 #undef __FUNCT__ 1587 #define __FUNCT__ "TSGetTimeStep" 1588 /*@ 1589 TSGetTimeStep - Gets the current timestep size. 1590 1591 Not Collective 1592 1593 Input Parameter: 1594 . ts - the TS context obtained from TSCreate() 1595 1596 Output Parameter: 1597 . dt - the current timestep size 1598 1599 Level: intermediate 1600 1601 .seealso: TSSetInitialTimeStep(), TSGetTimeStep() 1602 1603 .keywords: TS, get, timestep 1604 @*/ 1605 PetscErrorCode TSGetTimeStep(TS ts,PetscReal *dt) 1606 { 1607 PetscFunctionBegin; 1608 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1609 PetscValidRealPointer(dt,2); 1610 *dt = ts->time_step; 1611 PetscFunctionReturn(0); 1612 } 1613 1614 #undef __FUNCT__ 1615 #define __FUNCT__ "TSGetSolution" 1616 /*@ 1617 TSGetSolution - Returns the solution at the present timestep. It 1618 is valid to call this routine inside the function that you are evaluating 1619 in order to move to the new timestep. This vector not changed until 1620 the solution at the next timestep has been calculated. 1621 1622 Not Collective, but Vec returned is parallel if TS is parallel 1623 1624 Input Parameter: 1625 . ts - the TS context obtained from TSCreate() 1626 1627 Output Parameter: 1628 . v - the vector containing the solution 1629 1630 Level: intermediate 1631 1632 .seealso: TSGetTimeStep() 1633 1634 .keywords: TS, timestep, get, solution 1635 @*/ 1636 PetscErrorCode TSGetSolution(TS ts,Vec *v) 1637 { 1638 PetscFunctionBegin; 1639 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1640 PetscValidPointer(v,2); 1641 *v = ts->vec_sol; 1642 PetscFunctionReturn(0); 1643 } 1644 1645 /* ----- Routines to initialize and destroy a timestepper ---- */ 1646 #undef __FUNCT__ 1647 #define __FUNCT__ "TSSetProblemType" 1648 /*@ 1649 TSSetProblemType - Sets the type of problem to be solved. 1650 1651 Not collective 1652 1653 Input Parameters: 1654 + ts - The TS 1655 - type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1656 .vb 1657 U_t - A U = 0 (linear) 1658 U_t - A(t) U = 0 (linear) 1659 F(t,U,U_t) = 0 (nonlinear) 1660 .ve 1661 1662 Level: beginner 1663 1664 .keywords: TS, problem type 1665 .seealso: TSSetUp(), TSProblemType, TS 1666 @*/ 1667 PetscErrorCode TSSetProblemType(TS ts, TSProblemType type) 1668 { 1669 PetscErrorCode ierr; 1670 1671 PetscFunctionBegin; 1672 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1673 ts->problem_type = type; 1674 if (type == TS_LINEAR) { 1675 SNES snes; 1676 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1677 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 1678 } 1679 PetscFunctionReturn(0); 1680 } 1681 1682 #undef __FUNCT__ 1683 #define __FUNCT__ "TSGetProblemType" 1684 /*@C 1685 TSGetProblemType - Gets the type of problem to be solved. 1686 1687 Not collective 1688 1689 Input Parameter: 1690 . ts - The TS 1691 1692 Output Parameter: 1693 . type - One of TS_LINEAR, TS_NONLINEAR where these types refer to problems of the forms 1694 .vb 1695 M U_t = A U 1696 M(t) U_t = A(t) U 1697 F(t,U,U_t) 1698 .ve 1699 1700 Level: beginner 1701 1702 .keywords: TS, problem type 1703 .seealso: TSSetUp(), TSProblemType, TS 1704 @*/ 1705 PetscErrorCode TSGetProblemType(TS ts, TSProblemType *type) 1706 { 1707 PetscFunctionBegin; 1708 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 1709 PetscValidIntPointer(type,2); 1710 *type = ts->problem_type; 1711 PetscFunctionReturn(0); 1712 } 1713 1714 #undef __FUNCT__ 1715 #define __FUNCT__ "TSSetUp" 1716 /*@ 1717 TSSetUp - Sets up the internal data structures for the later use 1718 of a timestepper. 1719 1720 Collective on TS 1721 1722 Input Parameter: 1723 . ts - the TS context obtained from TSCreate() 1724 1725 Notes: 1726 For basic use of the TS solvers the user need not explicitly call 1727 TSSetUp(), since these actions will automatically occur during 1728 the call to TSStep(). However, if one wishes to control this 1729 phase separately, TSSetUp() should be called after TSCreate() 1730 and optional routines of the form TSSetXXX(), but before TSStep(). 1731 1732 Level: advanced 1733 1734 .keywords: TS, timestep, setup 1735 1736 .seealso: TSCreate(), TSStep(), TSDestroy() 1737 @*/ 1738 PetscErrorCode TSSetUp(TS ts) 1739 { 1740 PetscErrorCode ierr; 1741 DM dm; 1742 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 1743 PetscErrorCode (*jac)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 1744 TSIJacobian ijac; 1745 TSRHSJacobian rhsjac; 1746 1747 PetscFunctionBegin; 1748 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1749 if (ts->setupcalled) PetscFunctionReturn(0); 1750 1751 if (!((PetscObject)ts)->type_name) { 1752 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr); 1753 } 1754 1755 if (!ts->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call TSSetSolution() first"); 1756 1757 ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 1758 1759 if (ts->rhsjacobian.reuse) { 1760 Mat Amat,Pmat; 1761 SNES snes; 1762 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1763 ierr = SNESGetJacobian(snes,&Amat,&Pmat,NULL,NULL);CHKERRQ(ierr); 1764 /* Matching matrices implies that an IJacobian is NOT set, because if it had been set, the IJacobian's matrix would 1765 * have displaced the RHS matrix */ 1766 if (Amat == ts->Arhs) { 1767 ierr = MatDuplicate(ts->Arhs,MAT_DO_NOT_COPY_VALUES,&Amat);CHKERRQ(ierr); 1768 ierr = SNESSetJacobian(snes,Amat,NULL,NULL,NULL);CHKERRQ(ierr); 1769 ierr = MatDestroy(&Amat);CHKERRQ(ierr); 1770 } 1771 if (Pmat == ts->Brhs) { 1772 ierr = MatDuplicate(ts->Brhs,MAT_DO_NOT_COPY_VALUES,&Pmat);CHKERRQ(ierr); 1773 ierr = SNESSetJacobian(snes,NULL,Pmat,NULL,NULL);CHKERRQ(ierr); 1774 ierr = MatDestroy(&Pmat);CHKERRQ(ierr); 1775 } 1776 } 1777 1778 if (ts->ops->setup) { 1779 ierr = (*ts->ops->setup)(ts);CHKERRQ(ierr); 1780 } 1781 1782 /* in the case where we've set a DMTSFunction or what have you, we need the default SNESFunction 1783 to be set right but can't do it elsewhere due to the overreliance on ctx=ts. 1784 */ 1785 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 1786 ierr = DMSNESGetFunction(dm,&func,NULL);CHKERRQ(ierr); 1787 if (!func) { 1788 ierr =DMSNESSetFunction(dm,SNESTSFormFunction,ts);CHKERRQ(ierr); 1789 } 1790 /* if the SNES doesn't have a jacobian set and the TS has an ijacobian or rhsjacobian set, set the SNES to use it. 1791 Otherwise, the SNES will use coloring internally to form the Jacobian. 1792 */ 1793 ierr = DMSNESGetJacobian(dm,&jac,NULL);CHKERRQ(ierr); 1794 ierr = DMTSGetIJacobian(dm,&ijac,NULL);CHKERRQ(ierr); 1795 ierr = DMTSGetRHSJacobian(dm,&rhsjac,NULL);CHKERRQ(ierr); 1796 if (!jac && (ijac || rhsjac)) { 1797 ierr = DMSNESSetJacobian(dm,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1798 } 1799 ts->setupcalled = PETSC_TRUE; 1800 PetscFunctionReturn(0); 1801 } 1802 1803 #undef __FUNCT__ 1804 #define __FUNCT__ "TSReset" 1805 /*@ 1806 TSReset - Resets a TS context and removes any allocated Vecs and Mats. 1807 1808 Collective on TS 1809 1810 Input Parameter: 1811 . ts - the TS context obtained from TSCreate() 1812 1813 Level: beginner 1814 1815 .keywords: TS, timestep, reset 1816 1817 .seealso: TSCreate(), TSSetup(), TSDestroy() 1818 @*/ 1819 PetscErrorCode TSReset(TS ts) 1820 { 1821 PetscErrorCode ierr; 1822 1823 PetscFunctionBegin; 1824 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1825 if (ts->ops->reset) { 1826 ierr = (*ts->ops->reset)(ts);CHKERRQ(ierr); 1827 } 1828 if (ts->snes) {ierr = SNESReset(ts->snes);CHKERRQ(ierr);} 1829 1830 ierr = MatDestroy(&ts->Arhs);CHKERRQ(ierr); 1831 ierr = MatDestroy(&ts->Brhs);CHKERRQ(ierr); 1832 ierr = VecDestroy(&ts->Frhs);CHKERRQ(ierr); 1833 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 1834 ierr = VecDestroy(&ts->vatol);CHKERRQ(ierr); 1835 ierr = VecDestroy(&ts->vrtol);CHKERRQ(ierr); 1836 ierr = VecDestroyVecs(ts->nwork,&ts->work);CHKERRQ(ierr); 1837 1838 ts->setupcalled = PETSC_FALSE; 1839 PetscFunctionReturn(0); 1840 } 1841 1842 #undef __FUNCT__ 1843 #define __FUNCT__ "TSDestroy" 1844 /*@ 1845 TSDestroy - Destroys the timestepper context that was created 1846 with TSCreate(). 1847 1848 Collective on TS 1849 1850 Input Parameter: 1851 . ts - the TS context obtained from TSCreate() 1852 1853 Level: beginner 1854 1855 .keywords: TS, timestepper, destroy 1856 1857 .seealso: TSCreate(), TSSetUp(), TSSolve() 1858 @*/ 1859 PetscErrorCode TSDestroy(TS *ts) 1860 { 1861 PetscErrorCode ierr; 1862 1863 PetscFunctionBegin; 1864 if (!*ts) PetscFunctionReturn(0); 1865 PetscValidHeaderSpecific((*ts),TS_CLASSID,1); 1866 if (--((PetscObject)(*ts))->refct > 0) {*ts = 0; PetscFunctionReturn(0);} 1867 1868 ierr = TSReset((*ts));CHKERRQ(ierr); 1869 1870 /* if memory was published with AMS then destroy it */ 1871 ierr = PetscObjectAMSViewOff((PetscObject)*ts);CHKERRQ(ierr); 1872 if ((*ts)->ops->destroy) {ierr = (*(*ts)->ops->destroy)((*ts));CHKERRQ(ierr);} 1873 1874 ierr = TSAdaptDestroy(&(*ts)->adapt);CHKERRQ(ierr); 1875 ierr = SNESDestroy(&(*ts)->snes);CHKERRQ(ierr); 1876 ierr = DMDestroy(&(*ts)->dm);CHKERRQ(ierr); 1877 ierr = TSMonitorCancel((*ts));CHKERRQ(ierr); 1878 1879 ierr = PetscHeaderDestroy(ts);CHKERRQ(ierr); 1880 PetscFunctionReturn(0); 1881 } 1882 1883 #undef __FUNCT__ 1884 #define __FUNCT__ "TSGetSNES" 1885 /*@ 1886 TSGetSNES - Returns the SNES (nonlinear solver) associated with 1887 a TS (timestepper) context. Valid only for nonlinear problems. 1888 1889 Not Collective, but SNES is parallel if TS is parallel 1890 1891 Input Parameter: 1892 . ts - the TS context obtained from TSCreate() 1893 1894 Output Parameter: 1895 . snes - the nonlinear solver context 1896 1897 Notes: 1898 The user can then directly manipulate the SNES context to set various 1899 options, etc. Likewise, the user can then extract and manipulate the 1900 KSP, KSP, and PC contexts as well. 1901 1902 TSGetSNES() does not work for integrators that do not use SNES; in 1903 this case TSGetSNES() returns NULL in snes. 1904 1905 Level: beginner 1906 1907 .keywords: timestep, get, SNES 1908 @*/ 1909 PetscErrorCode TSGetSNES(TS ts,SNES *snes) 1910 { 1911 PetscErrorCode ierr; 1912 1913 PetscFunctionBegin; 1914 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1915 PetscValidPointer(snes,2); 1916 if (!ts->snes) { 1917 ierr = SNESCreate(PetscObjectComm((PetscObject)ts),&ts->snes);CHKERRQ(ierr); 1918 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 1919 ierr = PetscLogObjectParent(ts,ts->snes);CHKERRQ(ierr); 1920 ierr = PetscObjectIncrementTabLevel((PetscObject)ts->snes,(PetscObject)ts,1);CHKERRQ(ierr); 1921 if (ts->dm) {ierr = SNESSetDM(ts->snes,ts->dm);CHKERRQ(ierr);} 1922 if (ts->problem_type == TS_LINEAR) { 1923 ierr = SNESSetType(ts->snes,SNESKSPONLY);CHKERRQ(ierr); 1924 } 1925 } 1926 *snes = ts->snes; 1927 PetscFunctionReturn(0); 1928 } 1929 1930 #undef __FUNCT__ 1931 #define __FUNCT__ "TSSetSNES" 1932 /*@ 1933 TSSetSNES - Set the SNES (nonlinear solver) to be used by the timestepping context 1934 1935 Collective 1936 1937 Input Parameter: 1938 + ts - the TS context obtained from TSCreate() 1939 - snes - the nonlinear solver context 1940 1941 Notes: 1942 Most users should have the TS created by calling TSGetSNES() 1943 1944 Level: developer 1945 1946 .keywords: timestep, set, SNES 1947 @*/ 1948 PetscErrorCode TSSetSNES(TS ts,SNES snes) 1949 { 1950 PetscErrorCode ierr; 1951 PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 1952 1953 PetscFunctionBegin; 1954 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 1955 PetscValidHeaderSpecific(snes,SNES_CLASSID,2); 1956 ierr = PetscObjectReference((PetscObject)snes);CHKERRQ(ierr); 1957 ierr = SNESDestroy(&ts->snes);CHKERRQ(ierr); 1958 1959 ts->snes = snes; 1960 1961 ierr = SNESSetFunction(ts->snes,NULL,SNESTSFormFunction,ts);CHKERRQ(ierr); 1962 ierr = SNESGetJacobian(ts->snes,NULL,NULL,&func,NULL);CHKERRQ(ierr); 1963 if (func == SNESTSFormJacobian) { 1964 ierr = SNESSetJacobian(ts->snes,NULL,NULL,SNESTSFormJacobian,ts);CHKERRQ(ierr); 1965 } 1966 PetscFunctionReturn(0); 1967 } 1968 1969 #undef __FUNCT__ 1970 #define __FUNCT__ "TSGetKSP" 1971 /*@ 1972 TSGetKSP - Returns the KSP (linear solver) associated with 1973 a TS (timestepper) context. 1974 1975 Not Collective, but KSP is parallel if TS is parallel 1976 1977 Input Parameter: 1978 . ts - the TS context obtained from TSCreate() 1979 1980 Output Parameter: 1981 . ksp - the nonlinear solver context 1982 1983 Notes: 1984 The user can then directly manipulate the KSP context to set various 1985 options, etc. Likewise, the user can then extract and manipulate the 1986 KSP and PC contexts as well. 1987 1988 TSGetKSP() does not work for integrators that do not use KSP; 1989 in this case TSGetKSP() returns NULL in ksp. 1990 1991 Level: beginner 1992 1993 .keywords: timestep, get, KSP 1994 @*/ 1995 PetscErrorCode TSGetKSP(TS ts,KSP *ksp) 1996 { 1997 PetscErrorCode ierr; 1998 SNES snes; 1999 2000 PetscFunctionBegin; 2001 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2002 PetscValidPointer(ksp,2); 2003 if (!((PetscObject)ts)->type_name) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"KSP is not created yet. Call TSSetType() first"); 2004 if (ts->problem_type != TS_LINEAR) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Linear only; use TSGetSNES()"); 2005 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 2006 ierr = SNESGetKSP(snes,ksp);CHKERRQ(ierr); 2007 PetscFunctionReturn(0); 2008 } 2009 2010 /* ----------- Routines to set solver parameters ---------- */ 2011 2012 #undef __FUNCT__ 2013 #define __FUNCT__ "TSGetDuration" 2014 /*@ 2015 TSGetDuration - Gets the maximum number of timesteps to use and 2016 maximum time for iteration. 2017 2018 Not Collective 2019 2020 Input Parameters: 2021 + ts - the TS context obtained from TSCreate() 2022 . maxsteps - maximum number of iterations to use, or NULL 2023 - maxtime - final time to iterate to, or NULL 2024 2025 Level: intermediate 2026 2027 .keywords: TS, timestep, get, maximum, iterations, time 2028 @*/ 2029 PetscErrorCode TSGetDuration(TS ts, PetscInt *maxsteps, PetscReal *maxtime) 2030 { 2031 PetscFunctionBegin; 2032 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2033 if (maxsteps) { 2034 PetscValidIntPointer(maxsteps,2); 2035 *maxsteps = ts->max_steps; 2036 } 2037 if (maxtime) { 2038 PetscValidScalarPointer(maxtime,3); 2039 *maxtime = ts->max_time; 2040 } 2041 PetscFunctionReturn(0); 2042 } 2043 2044 #undef __FUNCT__ 2045 #define __FUNCT__ "TSSetDuration" 2046 /*@ 2047 TSSetDuration - Sets the maximum number of timesteps to use and 2048 maximum time for iteration. 2049 2050 Logically Collective on TS 2051 2052 Input Parameters: 2053 + ts - the TS context obtained from TSCreate() 2054 . maxsteps - maximum number of iterations to use 2055 - maxtime - final time to iterate to 2056 2057 Options Database Keys: 2058 . -ts_max_steps <maxsteps> - Sets maxsteps 2059 . -ts_final_time <maxtime> - Sets maxtime 2060 2061 Notes: 2062 The default maximum number of iterations is 5000. Default time is 5.0 2063 2064 Level: intermediate 2065 2066 .keywords: TS, timestep, set, maximum, iterations 2067 2068 .seealso: TSSetExactFinalTime() 2069 @*/ 2070 PetscErrorCode TSSetDuration(TS ts,PetscInt maxsteps,PetscReal maxtime) 2071 { 2072 PetscFunctionBegin; 2073 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2074 PetscValidLogicalCollectiveInt(ts,maxsteps,2); 2075 PetscValidLogicalCollectiveReal(ts,maxtime,2); 2076 if (maxsteps >= 0) ts->max_steps = maxsteps; 2077 if (maxtime != PETSC_DEFAULT) ts->max_time = maxtime; 2078 PetscFunctionReturn(0); 2079 } 2080 2081 #undef __FUNCT__ 2082 #define __FUNCT__ "TSSetSolution" 2083 /*@ 2084 TSSetSolution - Sets the initial solution vector 2085 for use by the TS routines. 2086 2087 Logically Collective on TS and Vec 2088 2089 Input Parameters: 2090 + ts - the TS context obtained from TSCreate() 2091 - u - the solution vector 2092 2093 Level: beginner 2094 2095 .keywords: TS, timestep, set, solution, initial conditions 2096 @*/ 2097 PetscErrorCode TSSetSolution(TS ts,Vec u) 2098 { 2099 PetscErrorCode ierr; 2100 DM dm; 2101 2102 PetscFunctionBegin; 2103 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2104 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2105 ierr = PetscObjectReference((PetscObject)u);CHKERRQ(ierr); 2106 ierr = VecDestroy(&ts->vec_sol);CHKERRQ(ierr); 2107 2108 ts->vec_sol = u; 2109 2110 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 2111 ierr = DMShellSetGlobalVector(dm,u);CHKERRQ(ierr); 2112 PetscFunctionReturn(0); 2113 } 2114 2115 #undef __FUNCT__ 2116 #define __FUNCT__ "TSSetPreStep" 2117 /*@C 2118 TSSetPreStep - Sets the general-purpose function 2119 called once at the beginning of each time step. 2120 2121 Logically Collective on TS 2122 2123 Input Parameters: 2124 + ts - The TS context obtained from TSCreate() 2125 - func - The function 2126 2127 Calling sequence of func: 2128 . func (TS ts); 2129 2130 Level: intermediate 2131 2132 Note: 2133 If a step is rejected, TSStep() will call this routine again before each attempt. 2134 The last completed time step number can be queried using TSGetTimeStepNumber(), the 2135 size of the step being attempted can be obtained using TSGetTimeStep(). 2136 2137 .keywords: TS, timestep 2138 .seealso: TSSetPreStage(), TSSetPostStep(), TSStep() 2139 @*/ 2140 PetscErrorCode TSSetPreStep(TS ts, PetscErrorCode (*func)(TS)) 2141 { 2142 PetscFunctionBegin; 2143 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2144 ts->prestep = func; 2145 PetscFunctionReturn(0); 2146 } 2147 2148 #undef __FUNCT__ 2149 #define __FUNCT__ "TSPreStep" 2150 /*@ 2151 TSPreStep - Runs the user-defined pre-step function. 2152 2153 Collective on TS 2154 2155 Input Parameters: 2156 . ts - The TS context obtained from TSCreate() 2157 2158 Notes: 2159 TSPreStep() is typically used within time stepping implementations, 2160 so most users would not generally call this routine themselves. 2161 2162 Level: developer 2163 2164 .keywords: TS, timestep 2165 .seealso: TSSetPreStep(), TSPreStage(), TSPostStep() 2166 @*/ 2167 PetscErrorCode TSPreStep(TS ts) 2168 { 2169 PetscErrorCode ierr; 2170 2171 PetscFunctionBegin; 2172 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2173 if (ts->prestep) { 2174 PetscStackCallStandard((*ts->prestep),(ts)); 2175 } 2176 PetscFunctionReturn(0); 2177 } 2178 2179 #undef __FUNCT__ 2180 #define __FUNCT__ "TSSetPreStage" 2181 /*@C 2182 TSSetPreStage - Sets the general-purpose function 2183 called once at the beginning of each stage. 2184 2185 Logically Collective on TS 2186 2187 Input Parameters: 2188 + ts - The TS context obtained from TSCreate() 2189 - func - The function 2190 2191 Calling sequence of func: 2192 . PetscErrorCode func(TS ts, PetscReal stagetime); 2193 2194 Level: intermediate 2195 2196 Note: 2197 There may be several stages per time step. If the solve for a given stage fails, the step may be rejected and retried. 2198 The time step number being computed can be queried using TSGetTimeStepNumber() and the total size of the step being 2199 attempted can be obtained using TSGetTimeStep(). The time at the start of the step is available via TSGetTime(). 2200 2201 .keywords: TS, timestep 2202 .seealso: TSSetPreStep(), TSSetPostStep(), TSGetApplicationContext() 2203 @*/ 2204 PetscErrorCode TSSetPreStage(TS ts, PetscErrorCode (*func)(TS,PetscReal)) 2205 { 2206 PetscFunctionBegin; 2207 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2208 ts->prestage = func; 2209 PetscFunctionReturn(0); 2210 } 2211 2212 #undef __FUNCT__ 2213 #define __FUNCT__ "TSPreStage" 2214 /*@ 2215 TSPreStage - Runs the user-defined pre-stage function set using TSSetPreStage() 2216 2217 Collective on TS 2218 2219 Input Parameters: 2220 . ts - The TS context obtained from TSCreate() 2221 2222 Notes: 2223 TSPreStage() is typically used within time stepping implementations, 2224 most users would not generally call this routine themselves. 2225 2226 Level: developer 2227 2228 .keywords: TS, timestep 2229 .seealso: TSSetPreStep(), TSPreStep(), TSPostStep() 2230 @*/ 2231 PetscErrorCode TSPreStage(TS ts, PetscReal stagetime) 2232 { 2233 PetscErrorCode ierr; 2234 2235 PetscFunctionBegin; 2236 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2237 if (ts->prestage) { 2238 PetscStackCallStandard((*ts->prestage),(ts,stagetime)); 2239 } 2240 PetscFunctionReturn(0); 2241 } 2242 2243 #undef __FUNCT__ 2244 #define __FUNCT__ "TSSetPostStep" 2245 /*@C 2246 TSSetPostStep - Sets the general-purpose function 2247 called once at the end of each time step. 2248 2249 Logically Collective on TS 2250 2251 Input Parameters: 2252 + ts - The TS context obtained from TSCreate() 2253 - func - The function 2254 2255 Calling sequence of func: 2256 $ func (TS ts); 2257 2258 Level: intermediate 2259 2260 .keywords: TS, timestep 2261 .seealso: TSSetPreStep(), TSSetPreStage(), TSGetTimeStep(), TSGetTimeStepNumber(), TSGetTime() 2262 @*/ 2263 PetscErrorCode TSSetPostStep(TS ts, PetscErrorCode (*func)(TS)) 2264 { 2265 PetscFunctionBegin; 2266 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2267 ts->poststep = func; 2268 PetscFunctionReturn(0); 2269 } 2270 2271 #undef __FUNCT__ 2272 #define __FUNCT__ "TSPostStep" 2273 /*@ 2274 TSPostStep - Runs the user-defined post-step function. 2275 2276 Collective on TS 2277 2278 Input Parameters: 2279 . ts - The TS context obtained from TSCreate() 2280 2281 Notes: 2282 TSPostStep() is typically used within time stepping implementations, 2283 so most users would not generally call this routine themselves. 2284 2285 Level: developer 2286 2287 .keywords: TS, timestep 2288 @*/ 2289 PetscErrorCode TSPostStep(TS ts) 2290 { 2291 PetscErrorCode ierr; 2292 2293 PetscFunctionBegin; 2294 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2295 if (ts->poststep) { 2296 PetscStackCallStandard((*ts->poststep),(ts)); 2297 } 2298 PetscFunctionReturn(0); 2299 } 2300 2301 /* ------------ Routines to set performance monitoring options ----------- */ 2302 2303 #undef __FUNCT__ 2304 #define __FUNCT__ "TSMonitorSet" 2305 /*@C 2306 TSMonitorSet - Sets an ADDITIONAL function that is to be used at every 2307 timestep to display the iteration's progress. 2308 2309 Logically Collective on TS 2310 2311 Input Parameters: 2312 + ts - the TS context obtained from TSCreate() 2313 . monitor - monitoring routine 2314 . mctx - [optional] user-defined context for private data for the 2315 monitor routine (use NULL if no context is desired) 2316 - monitordestroy - [optional] routine that frees monitor context 2317 (may be NULL) 2318 2319 Calling sequence of monitor: 2320 $ int monitor(TS ts,PetscInt steps,PetscReal time,Vec u,void *mctx) 2321 2322 + ts - the TS context 2323 . 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 2324 been interpolated to) 2325 . time - current time 2326 . u - current iterate 2327 - mctx - [optional] monitoring context 2328 2329 Notes: 2330 This routine adds an additional monitor to the list of monitors that 2331 already has been loaded. 2332 2333 Fortran notes: Only a single monitor function can be set for each TS object 2334 2335 Level: intermediate 2336 2337 .keywords: TS, timestep, set, monitor 2338 2339 .seealso: TSMonitorDefault(), TSMonitorCancel() 2340 @*/ 2341 PetscErrorCode TSMonitorSet(TS ts,PetscErrorCode (*monitor)(TS,PetscInt,PetscReal,Vec,void*),void *mctx,PetscErrorCode (*mdestroy)(void**)) 2342 { 2343 PetscFunctionBegin; 2344 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2345 if (ts->numbermonitors >= MAXTSMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 2346 ts->monitor[ts->numbermonitors] = monitor; 2347 ts->monitordestroy[ts->numbermonitors] = mdestroy; 2348 ts->monitorcontext[ts->numbermonitors++] = (void*)mctx; 2349 PetscFunctionReturn(0); 2350 } 2351 2352 #undef __FUNCT__ 2353 #define __FUNCT__ "TSMonitorCancel" 2354 /*@C 2355 TSMonitorCancel - Clears all the monitors that have been set on a time-step object. 2356 2357 Logically Collective on TS 2358 2359 Input Parameters: 2360 . ts - the TS context obtained from TSCreate() 2361 2362 Notes: 2363 There is no way to remove a single, specific monitor. 2364 2365 Level: intermediate 2366 2367 .keywords: TS, timestep, set, monitor 2368 2369 .seealso: TSMonitorDefault(), TSMonitorSet() 2370 @*/ 2371 PetscErrorCode TSMonitorCancel(TS ts) 2372 { 2373 PetscErrorCode ierr; 2374 PetscInt i; 2375 2376 PetscFunctionBegin; 2377 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2378 for (i=0; i<ts->numbermonitors; i++) { 2379 if (ts->monitordestroy[i]) { 2380 ierr = (*ts->monitordestroy[i])(&ts->monitorcontext[i]);CHKERRQ(ierr); 2381 } 2382 } 2383 ts->numbermonitors = 0; 2384 PetscFunctionReturn(0); 2385 } 2386 2387 #undef __FUNCT__ 2388 #define __FUNCT__ "TSMonitorDefault" 2389 /*@ 2390 TSMonitorDefault - Sets the Default monitor 2391 2392 Level: intermediate 2393 2394 .keywords: TS, set, monitor 2395 2396 .seealso: TSMonitorDefault(), TSMonitorSet() 2397 @*/ 2398 PetscErrorCode TSMonitorDefault(TS ts,PetscInt step,PetscReal ptime,Vec v,void *dummy) 2399 { 2400 PetscErrorCode ierr; 2401 PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ts)); 2402 2403 PetscFunctionBegin; 2404 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2405 ierr = PetscViewerASCIIPrintf(viewer,"%D TS dt %g time %g\n",step,(double)ts->time_step,(double)ptime);CHKERRQ(ierr); 2406 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ts)->tablevel);CHKERRQ(ierr); 2407 PetscFunctionReturn(0); 2408 } 2409 2410 #undef __FUNCT__ 2411 #define __FUNCT__ "TSSetRetainStages" 2412 /*@ 2413 TSSetRetainStages - Request that all stages in the upcoming step be stored so that interpolation will be available. 2414 2415 Logically Collective on TS 2416 2417 Input Argument: 2418 . ts - time stepping context 2419 2420 Output Argument: 2421 . flg - PETSC_TRUE or PETSC_FALSE 2422 2423 Level: intermediate 2424 2425 .keywords: TS, set 2426 2427 .seealso: TSInterpolate(), TSSetPostStep() 2428 @*/ 2429 PetscErrorCode TSSetRetainStages(TS ts,PetscBool flg) 2430 { 2431 PetscFunctionBegin; 2432 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2433 ts->retain_stages = flg; 2434 PetscFunctionReturn(0); 2435 } 2436 2437 #undef __FUNCT__ 2438 #define __FUNCT__ "TSInterpolate" 2439 /*@ 2440 TSInterpolate - Interpolate the solution computed during the previous step to an arbitrary location in the interval 2441 2442 Collective on TS 2443 2444 Input Argument: 2445 + ts - time stepping context 2446 - t - time to interpolate to 2447 2448 Output Argument: 2449 . U - state at given time 2450 2451 Notes: 2452 The user should call TSSetRetainStages() before taking a step in which interpolation will be requested. 2453 2454 Level: intermediate 2455 2456 Developer Notes: 2457 TSInterpolate() and the storing of previous steps/stages should be generalized to support delay differential equations and continuous adjoints. 2458 2459 .keywords: TS, set 2460 2461 .seealso: TSSetRetainStages(), TSSetPostStep() 2462 @*/ 2463 PetscErrorCode TSInterpolate(TS ts,PetscReal t,Vec U) 2464 { 2465 PetscErrorCode ierr; 2466 2467 PetscFunctionBegin; 2468 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2469 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 2470 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); 2471 if (!ts->ops->interpolate) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"%s does not provide interpolation",((PetscObject)ts)->type_name); 2472 ierr = (*ts->ops->interpolate)(ts,t,U);CHKERRQ(ierr); 2473 PetscFunctionReturn(0); 2474 } 2475 2476 #undef __FUNCT__ 2477 #define __FUNCT__ "TSStep" 2478 /*@ 2479 TSStep - Steps one time step 2480 2481 Collective on TS 2482 2483 Input Parameter: 2484 . ts - the TS context obtained from TSCreate() 2485 2486 Level: intermediate 2487 2488 Notes: 2489 The hook set using TSSetPreStep() is called before each attempt to take the step. In general, the time step size may 2490 be changed due to adaptive error controller or solve failures. Note that steps may contain multiple stages. 2491 2492 This may over-step the final time provided in TSSetDuration() depending on the time-step used. TSSolve() interpolates to exactly the 2493 time provided in TSSetDuration(). One can use TSInterpolate() to determine an interpolated solution within the final timestep. 2494 2495 .keywords: TS, timestep, solve 2496 2497 .seealso: TSCreate(), TSSetUp(), TSDestroy(), TSSolve(), TSSetPreStep(), TSSetPreStage(), TSInterpolate() 2498 @*/ 2499 PetscErrorCode TSStep(TS ts) 2500 { 2501 PetscReal ptime_prev; 2502 PetscErrorCode ierr; 2503 2504 PetscFunctionBegin; 2505 PetscValidHeaderSpecific(ts, TS_CLASSID,1); 2506 ierr = TSSetUp(ts);CHKERRQ(ierr); 2507 2508 ts->reason = TS_CONVERGED_ITERATING; 2509 ptime_prev = ts->ptime; 2510 2511 ierr = PetscLogEventBegin(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2512 ierr = (*ts->ops->step)(ts);CHKERRQ(ierr); 2513 ierr = PetscLogEventEnd(TS_Step,ts,0,0,0);CHKERRQ(ierr); 2514 2515 ts->time_step_prev = ts->ptime - ptime_prev; 2516 2517 if (ts->reason < 0) { 2518 if (ts->errorifstepfailed) { 2519 if (ts->reason == TS_DIVERGED_NONLINEAR_SOLVE) { 2520 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]); 2521 } else SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_NOT_CONVERGED,"TSStep has failed due to %s",TSConvergedReasons[ts->reason]); 2522 } 2523 } else if (!ts->reason) { 2524 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 2525 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 2526 } 2527 PetscFunctionReturn(0); 2528 } 2529 2530 #undef __FUNCT__ 2531 #define __FUNCT__ "TSEvaluateStep" 2532 /*@ 2533 TSEvaluateStep - Evaluate the solution at the end of a time step with a given order of accuracy. 2534 2535 Collective on TS 2536 2537 Input Arguments: 2538 + ts - time stepping context 2539 . order - desired order of accuracy 2540 - done - whether the step was evaluated at this order (pass NULL to generate an error if not available) 2541 2542 Output Arguments: 2543 . U - state at the end of the current step 2544 2545 Level: advanced 2546 2547 Notes: 2548 This function cannot be called until all stages have been evaluated. 2549 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. 2550 2551 .seealso: TSStep(), TSAdapt 2552 @*/ 2553 PetscErrorCode TSEvaluateStep(TS ts,PetscInt order,Vec U,PetscBool *done) 2554 { 2555 PetscErrorCode ierr; 2556 2557 PetscFunctionBegin; 2558 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2559 PetscValidType(ts,1); 2560 PetscValidHeaderSpecific(U,VEC_CLASSID,3); 2561 if (!ts->ops->evaluatestep) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_SUP,"TSEvaluateStep not implemented for type '%s'",((PetscObject)ts)->type_name); 2562 ierr = (*ts->ops->evaluatestep)(ts,order,U,done);CHKERRQ(ierr); 2563 PetscFunctionReturn(0); 2564 } 2565 2566 #undef __FUNCT__ 2567 #define __FUNCT__ "TSSolve" 2568 /*@ 2569 TSSolve - Steps the requested number of timesteps. 2570 2571 Collective on TS 2572 2573 Input Parameter: 2574 + ts - the TS context obtained from TSCreate() 2575 - u - the solution vector (can be null if TSSetSolution() was used, otherwise must contain the initial conditions) 2576 2577 Level: beginner 2578 2579 Notes: 2580 The final time returned by this function may be different from the time of the internally 2581 held state accessible by TSGetSolution() and TSGetTime() because the method may have 2582 stepped over the final time. 2583 2584 .keywords: TS, timestep, solve 2585 2586 .seealso: TSCreate(), TSSetSolution(), TSStep() 2587 @*/ 2588 PetscErrorCode TSSolve(TS ts,Vec u) 2589 { 2590 PetscBool flg; 2591 PetscViewer viewer; 2592 Vec solution; 2593 PetscErrorCode ierr; 2594 PetscViewerFormat format; 2595 2596 PetscFunctionBegin; 2597 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 2598 if (u) PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2599 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 */ 2600 PetscValidHeaderSpecific(u,VEC_CLASSID,2); 2601 if (!ts->vec_sol || u == ts->vec_sol) { 2602 ierr = VecDuplicate(u,&solution);CHKERRQ(ierr); 2603 ierr = TSSetSolution(ts,solution);CHKERRQ(ierr); 2604 ierr = VecDestroy(&solution);CHKERRQ(ierr); /* grant ownership */ 2605 } 2606 ierr = VecCopy(u,ts->vec_sol);CHKERRQ(ierr); 2607 } else if (u) { 2608 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 2609 } 2610 ierr = TSSetUp(ts);CHKERRQ(ierr); 2611 /* reset time step and iteration counters */ 2612 ts->steps = 0; 2613 ts->ksp_its = 0; 2614 ts->snes_its = 0; 2615 ts->num_snes_failures = 0; 2616 ts->reject = 0; 2617 ts->reason = TS_CONVERGED_ITERATING; 2618 2619 ierr = PetscOptionsGetViewer(PetscObjectComm((PetscObject)ts),((PetscObject)ts)->prefix,"-ts_view_pre",&viewer,&format,&flg);CHKERRQ(ierr); 2620 if (flg && !PetscPreLoadingOn) { 2621 ierr = PetscViewerPushFormat(viewer,format);CHKERRQ(ierr); 2622 ierr = TSView(ts,viewer);CHKERRQ(ierr); 2623 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 2624 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 2625 } 2626 2627 if (ts->ops->solve) { /* This private interface is transitional and should be removed when all implementations are updated. */ 2628 ierr = (*ts->ops->solve)(ts);CHKERRQ(ierr); 2629 ierr = VecCopy(ts->vec_sol,u);CHKERRQ(ierr); 2630 ts->solvetime = ts->ptime; 2631 } else { 2632 /* steps the requested number of timesteps. */ 2633 if (ts->steps >= ts->max_steps) ts->reason = TS_CONVERGED_ITS; 2634 else if (ts->ptime >= ts->max_time) ts->reason = TS_CONVERGED_TIME; 2635 while (!ts->reason) { 2636 ierr = TSMonitor(ts,ts->steps,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 2637 ierr = TSStep(ts);CHKERRQ(ierr); 2638 ierr = TSPostStep(ts);CHKERRQ(ierr); 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