1 /* 2 Code for timestepping with BDF methods 3 */ 4 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 5 #include <petscdm.h> 6 7 static PetscBool cited = PETSC_FALSE; 8 static const char citation[] = 9 "@book{Brenan1995,\n" 10 " title = {Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations},\n" 11 " author = {Brenan, K. and Campbell, S. and Petzold, L.},\n" 12 " publisher = {Society for Industrial and Applied Mathematics},\n" 13 " year = {1995},\n" 14 " doi = {10.1137/1.9781611971224},\n}\n"; 15 16 typedef struct { 17 PetscInt k,n; 18 PetscReal time[6+2]; 19 Vec work[6+2]; 20 Vec tvwork[6+2]; 21 PetscReal shift; 22 Vec vec_dot; /* Xdot when !transientvar, else Cdot where C(X) is the transient variable. */ 23 Vec vec_wrk; 24 Vec vec_lte; 25 26 PetscBool transientvar; 27 PetscInt order; 28 TSStepStatus status; 29 } TS_BDF; 30 31 /* Compute Lagrange polynomials on T[:n] evaluated at t. 32 * If one has data (T[i], Y[i]), then the interpolation/extrapolation is f(t) = \sum_i L[i]*Y[i]. 33 */ 34 static inline void LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[]) 35 { 36 PetscInt k,j; 37 for (k=0; k<n; k++) 38 for (L[k]=1, j=0; j<n; j++) 39 if (j != k) 40 L[k] *= (t - T[j])/(T[k] - T[j]); 41 } 42 43 static inline void LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[]) 44 { 45 PetscInt k,j,i; 46 for (k=0; k<n; k++) 47 for (dL[k]=0, j=0; j<n; j++) 48 if (j != k) { 49 PetscReal L = 1/(T[k] - T[j]); 50 for (i=0; i<n; i++) 51 if (i != j && i != k) 52 L *= (t - T[i])/(T[k] - T[i]); 53 dL[k] += L; 54 } 55 } 56 57 static PetscErrorCode TSBDF_GetVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot) 58 { 59 TS_BDF *bdf = (TS_BDF*)ts->data; 60 61 PetscFunctionBegin; 62 if (dm && dm != ts->dm) { 63 CHKERRQ(DMGetNamedGlobalVector(dm,"TSBDF_Vec_Xdot",Xdot)); 64 CHKERRQ(DMGetNamedGlobalVector(dm,"TSBDF_Vec_Ydot",Ydot)); 65 } else { 66 *Xdot = bdf->vec_dot; 67 *Ydot = bdf->vec_wrk; 68 } 69 PetscFunctionReturn(0); 70 } 71 72 static PetscErrorCode TSBDF_RestoreVecs(TS ts,DM dm,Vec *Xdot,Vec *Ydot) 73 { 74 TS_BDF *bdf = (TS_BDF*)ts->data; 75 76 PetscFunctionBegin; 77 if (dm && dm != ts->dm) { 78 CHKERRQ(DMRestoreNamedGlobalVector(dm,"TSBDF_Vec_Xdot",Xdot)); 79 CHKERRQ(DMRestoreNamedGlobalVector(dm,"TSBDF_Vec_Ydot",Ydot)); 80 } else { 81 PetscCheck(*Xdot == bdf->vec_dot,PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Vec does not match the cache"); 82 PetscCheck(*Ydot == bdf->vec_wrk,PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_INCOMP,"Vec does not match the cache"); 83 *Xdot = NULL; 84 *Ydot = NULL; 85 } 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode DMCoarsenHook_TSBDF(DM fine,DM coarse,void *ctx) 90 { 91 PetscFunctionBegin; 92 PetscFunctionReturn(0); 93 } 94 95 static PetscErrorCode DMRestrictHook_TSBDF(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 96 { 97 TS ts = (TS)ctx; 98 Vec Ydot,Ydot_c; 99 Vec Xdot,Xdot_c; 100 101 PetscFunctionBegin; 102 CHKERRQ(TSBDF_GetVecs(ts,fine,&Xdot,&Ydot)); 103 CHKERRQ(TSBDF_GetVecs(ts,coarse,&Xdot_c,&Ydot_c)); 104 105 CHKERRQ(MatRestrict(restrct,Ydot,Ydot_c)); 106 CHKERRQ(VecPointwiseMult(Ydot_c,rscale,Ydot_c)); 107 108 CHKERRQ(TSBDF_RestoreVecs(ts,fine,&Xdot,&Ydot)); 109 CHKERRQ(TSBDF_RestoreVecs(ts,coarse,&Xdot_c,&Ydot_c)); 110 PetscFunctionReturn(0); 111 } 112 113 static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X) 114 { 115 TS_BDF *bdf = (TS_BDF*)ts->data; 116 PetscInt i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec)); 117 Vec tail = bdf->work[n-1],tvtail = bdf->tvwork[n-1]; 118 119 PetscFunctionBegin; 120 for (i=n-1; i>=2; i--) { 121 bdf->time[i] = bdf->time[i-1]; 122 bdf->work[i] = bdf->work[i-1]; 123 bdf->tvwork[i] = bdf->tvwork[i-1]; 124 } 125 bdf->n = PetscMin(bdf->n+1,n-1); 126 bdf->time[1] = t; 127 bdf->work[1] = tail; 128 bdf->tvwork[1] = tvtail; 129 CHKERRQ(VecCopy(X,tail)); 130 CHKERRQ(TSComputeTransientVariable(ts,tail,tvtail)); 131 PetscFunctionReturn(0); 132 } 133 134 static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte) 135 { 136 TS_BDF *bdf = (TS_BDF*)ts->data; 137 PetscInt i,n = order+1; 138 PetscReal *time = bdf->time; 139 Vec *vecs = bdf->work; 140 PetscScalar a[8],b[8],alpha[8]; 141 142 PetscFunctionBegin; 143 LagrangeBasisDers(n+0,time[0],time,a); a[n] =0; 144 LagrangeBasisDers(n+1,time[0],time,b); 145 for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0]; 146 CHKERRQ(VecZeroEntries(lte)); 147 CHKERRQ(VecMAXPY(lte,n+1,alpha,vecs)); 148 PetscFunctionReturn(0); 149 } 150 151 static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X) 152 { 153 TS_BDF *bdf = (TS_BDF*)ts->data; 154 PetscInt n = order+1; 155 PetscReal *time = bdf->time+1; 156 Vec *vecs = bdf->work+1; 157 PetscScalar alpha[7]; 158 159 PetscFunctionBegin; 160 n = PetscMin(n,bdf->n); 161 LagrangeBasisVals(n,t,time,alpha); 162 CHKERRQ(VecZeroEntries(X)); 163 CHKERRQ(VecMAXPY(X,n,alpha,vecs)); 164 PetscFunctionReturn(0); 165 } 166 167 static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X) 168 { 169 TS_BDF *bdf = (TS_BDF*)ts->data; 170 PetscInt n = order+1; 171 PetscReal *time = bdf->time; 172 Vec *vecs = bdf->work; 173 PetscScalar alpha[7]; 174 175 PetscFunctionBegin; 176 LagrangeBasisVals(n,t,time,alpha); 177 CHKERRQ(VecZeroEntries(X)); 178 CHKERRQ(VecMAXPY(X,n,alpha,vecs)); 179 PetscFunctionReturn(0); 180 } 181 182 /* Compute the affine term V0 such that Xdot = shift*X + V0. 183 * 184 * When using transient variables, we're computing Cdot = shift*C(X) + V0, and thus choose a linear combination of tvwork. 185 */ 186 static PetscErrorCode TSBDF_PreSolve(TS ts) 187 { 188 TS_BDF *bdf = (TS_BDF*)ts->data; 189 PetscInt i,n = PetscMax(bdf->k,1) + 1; 190 Vec V,V0; 191 Vec vecs[7]; 192 PetscScalar alpha[7]; 193 194 PetscFunctionBegin; 195 CHKERRQ(TSBDF_GetVecs(ts,NULL,&V,&V0)); 196 LagrangeBasisDers(n,bdf->time[0],bdf->time,alpha); 197 for (i=1; i<n; i++) { 198 vecs[i] = bdf->transientvar ? bdf->tvwork[i] : bdf->work[i]; 199 } 200 CHKERRQ(VecZeroEntries(V0)); 201 CHKERRQ(VecMAXPY(V0,n-1,alpha+1,vecs+1)); 202 bdf->shift = PetscRealPart(alpha[0]); 203 CHKERRQ(TSBDF_RestoreVecs(ts,NULL,&V,&V0)); 204 PetscFunctionReturn(0); 205 } 206 207 static PetscErrorCode TSBDF_SNESSolve(TS ts,Vec b,Vec x) 208 { 209 PetscInt nits,lits; 210 211 PetscFunctionBegin; 212 CHKERRQ(TSBDF_PreSolve(ts)); 213 CHKERRQ(SNESSolve(ts->snes,b,x)); 214 CHKERRQ(SNESGetIterationNumber(ts->snes,&nits)); 215 CHKERRQ(SNESGetLinearSolveIterations(ts->snes,&lits)); 216 ts->snes_its += nits; ts->ksp_its += lits; 217 PetscFunctionReturn(0); 218 } 219 220 static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept) 221 { 222 TS_BDF *bdf = (TS_BDF*)ts->data; 223 224 PetscFunctionBegin; 225 bdf->k = 1; bdf->n = 0; 226 CHKERRQ(TSBDF_Advance(ts,ts->ptime,ts->vec_sol)); 227 228 bdf->time[0] = ts->ptime + ts->time_step/2; 229 CHKERRQ(VecCopy(bdf->work[1],bdf->work[0])); 230 CHKERRQ(TSPreStage(ts,bdf->time[0])); 231 CHKERRQ(TSBDF_SNESSolve(ts,NULL,bdf->work[0])); 232 CHKERRQ(TSPostStage(ts,bdf->time[0],0,&bdf->work[0])); 233 CHKERRQ(TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept)); 234 if (!*accept) PetscFunctionReturn(0); 235 236 bdf->k = PetscMin(2,bdf->order); bdf->n++; 237 CHKERRQ(VecCopy(bdf->work[0],bdf->work[2])); 238 bdf->time[2] = bdf->time[0]; 239 CHKERRQ(TSComputeTransientVariable(ts,bdf->work[2],bdf->tvwork[2])); 240 PetscFunctionReturn(0); 241 } 242 243 static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"}; 244 245 static PetscErrorCode TSStep_BDF(TS ts) 246 { 247 TS_BDF *bdf = (TS_BDF*)ts->data; 248 PetscInt rejections = 0; 249 PetscBool stageok,accept = PETSC_TRUE; 250 PetscReal next_time_step = ts->time_step; 251 252 PetscFunctionBegin; 253 CHKERRQ(PetscCitationsRegister(citation,&cited)); 254 255 if (!ts->steprollback && !ts->steprestart) { 256 bdf->k = PetscMin(bdf->k+1,bdf->order); 257 CHKERRQ(TSBDF_Advance(ts,ts->ptime,ts->vec_sol)); 258 } 259 260 bdf->status = TS_STEP_INCOMPLETE; 261 while (!ts->reason && bdf->status != TS_STEP_COMPLETE) { 262 263 if (ts->steprestart) { 264 CHKERRQ(TSBDF_Restart(ts,&stageok)); 265 if (!stageok) goto reject_step; 266 } 267 268 bdf->time[0] = ts->ptime + ts->time_step; 269 CHKERRQ(TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0])); 270 CHKERRQ(TSPreStage(ts,bdf->time[0])); 271 CHKERRQ(TSBDF_SNESSolve(ts,NULL,bdf->work[0])); 272 CHKERRQ(TSPostStage(ts,bdf->time[0],0,&bdf->work[0])); 273 CHKERRQ(TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok)); 274 if (!stageok) goto reject_step; 275 276 bdf->status = TS_STEP_PENDING; 277 CHKERRQ(TSAdaptCandidatesClear(ts->adapt)); 278 CHKERRQ(TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE)); 279 CHKERRQ(TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept)); 280 bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 281 if (!accept) { ts->time_step = next_time_step; goto reject_step; } 282 283 CHKERRQ(VecCopy(bdf->work[0],ts->vec_sol)); 284 ts->ptime += ts->time_step; 285 ts->time_step = next_time_step; 286 break; 287 288 reject_step: 289 ts->reject++; accept = PETSC_FALSE; 290 if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 291 CHKERRQ(PetscInfo(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections)); 292 ts->reason = TS_DIVERGED_STEP_REJECTED; 293 } 294 } 295 PetscFunctionReturn(0); 296 } 297 298 static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X) 299 { 300 TS_BDF *bdf = (TS_BDF*)ts->data; 301 302 PetscFunctionBegin; 303 CHKERRQ(TSBDF_Interpolate(ts,bdf->k,t,X)); 304 PetscFunctionReturn(0); 305 } 306 307 static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte) 308 { 309 TS_BDF *bdf = (TS_BDF*)ts->data; 310 PetscInt k = bdf->k; 311 PetscReal wltea,wlter; 312 Vec X = bdf->work[0], Y = bdf->vec_lte; 313 314 PetscFunctionBegin; 315 k = PetscMin(k,bdf->n-1); 316 CHKERRQ(TSBDF_VecLTE(ts,k,Y)); 317 CHKERRQ(VecAXPY(Y,1,X)); 318 CHKERRQ(TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter)); 319 if (order) *order = k + 1; 320 PetscFunctionReturn(0); 321 } 322 323 static PetscErrorCode TSRollBack_BDF(TS ts) 324 { 325 TS_BDF *bdf = (TS_BDF*)ts->data; 326 327 PetscFunctionBegin; 328 CHKERRQ(VecCopy(bdf->work[1],ts->vec_sol)); 329 PetscFunctionReturn(0); 330 } 331 332 static PetscErrorCode SNESTSFormFunction_BDF(SNES snes,Vec X,Vec F,TS ts) 333 { 334 TS_BDF *bdf = (TS_BDF*)ts->data; 335 DM dm, dmsave = ts->dm; 336 PetscReal t = bdf->time[0]; 337 PetscReal shift = bdf->shift; 338 Vec V,V0; 339 340 PetscFunctionBegin; 341 CHKERRQ(SNESGetDM(snes,&dm)); 342 CHKERRQ(TSBDF_GetVecs(ts,dm,&V,&V0)); 343 if (bdf->transientvar) { /* shift*C(X) + V0 */ 344 CHKERRQ(TSComputeTransientVariable(ts,X,V)); 345 CHKERRQ(VecAYPX(V,shift,V0)); 346 } else { /* shift*X + V0 */ 347 CHKERRQ(VecWAXPY(V,shift,X,V0)); 348 } 349 350 /* F = Function(t,X,V) */ 351 ts->dm = dm; 352 CHKERRQ(TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE)); 353 ts->dm = dmsave; 354 355 CHKERRQ(TSBDF_RestoreVecs(ts,dm,&V,&V0)); 356 PetscFunctionReturn(0); 357 } 358 359 static PetscErrorCode SNESTSFormJacobian_BDF(SNES snes,Vec X,Mat J,Mat P,TS ts) 360 { 361 TS_BDF *bdf = (TS_BDF*)ts->data; 362 DM dm, dmsave = ts->dm; 363 PetscReal t = bdf->time[0]; 364 PetscReal shift = bdf->shift; 365 Vec V,V0; 366 367 PetscFunctionBegin; 368 CHKERRQ(SNESGetDM(snes,&dm)); 369 CHKERRQ(TSBDF_GetVecs(ts,dm,&V,&V0)); 370 371 /* J,P = Jacobian(t,X,V) */ 372 ts->dm = dm; 373 CHKERRQ(TSComputeIJacobian(ts,t,X,V,shift,J,P,PETSC_FALSE)); 374 ts->dm = dmsave; 375 376 CHKERRQ(TSBDF_RestoreVecs(ts,dm,&V,&V0)); 377 PetscFunctionReturn(0); 378 } 379 380 static PetscErrorCode TSReset_BDF(TS ts) 381 { 382 TS_BDF *bdf = (TS_BDF*)ts->data; 383 size_t i,n = sizeof(bdf->work)/sizeof(Vec); 384 385 PetscFunctionBegin; 386 bdf->k = bdf->n = 0; 387 for (i=0; i<n; i++) { 388 CHKERRQ(VecDestroy(&bdf->work[i])); 389 CHKERRQ(VecDestroy(&bdf->tvwork[i])); 390 } 391 CHKERRQ(VecDestroy(&bdf->vec_dot)); 392 CHKERRQ(VecDestroy(&bdf->vec_wrk)); 393 CHKERRQ(VecDestroy(&bdf->vec_lte)); 394 if (ts->dm) CHKERRQ(DMCoarsenHookRemove(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts)); 395 PetscFunctionReturn(0); 396 } 397 398 static PetscErrorCode TSDestroy_BDF(TS ts) 399 { 400 PetscFunctionBegin; 401 CHKERRQ(TSReset_BDF(ts)); 402 CHKERRQ(PetscFree(ts->data)); 403 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL)); 404 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL)); 405 PetscFunctionReturn(0); 406 } 407 408 static PetscErrorCode TSSetUp_BDF(TS ts) 409 { 410 TS_BDF *bdf = (TS_BDF*)ts->data; 411 size_t i,n = sizeof(bdf->work)/sizeof(Vec); 412 PetscReal low,high,two = 2; 413 414 PetscFunctionBegin; 415 CHKERRQ(TSHasTransientVariable(ts,&bdf->transientvar)); 416 bdf->k = bdf->n = 0; 417 for (i=0; i<n; i++) { 418 CHKERRQ(VecDuplicate(ts->vec_sol,&bdf->work[i])); 419 if (i && bdf->transientvar) { 420 CHKERRQ(VecDuplicate(ts->vec_sol,&bdf->tvwork[i])); 421 } 422 } 423 CHKERRQ(VecDuplicate(ts->vec_sol,&bdf->vec_dot)); 424 CHKERRQ(VecDuplicate(ts->vec_sol,&bdf->vec_wrk)); 425 CHKERRQ(VecDuplicate(ts->vec_sol,&bdf->vec_lte)); 426 CHKERRQ(TSGetDM(ts,&ts->dm)); 427 CHKERRQ(DMCoarsenHookAdd(ts->dm,DMCoarsenHook_TSBDF,DMRestrictHook_TSBDF,ts)); 428 429 CHKERRQ(TSGetAdapt(ts,&ts->adapt)); 430 CHKERRQ(TSAdaptCandidatesClear(ts->adapt)); 431 CHKERRQ(TSAdaptGetClip(ts->adapt,&low,&high)); 432 CHKERRQ(TSAdaptSetClip(ts->adapt,low,PetscMin(high,two))); 433 434 CHKERRQ(TSGetSNES(ts,&ts->snes)); 435 PetscFunctionReturn(0); 436 } 437 438 static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts) 439 { 440 PetscFunctionBegin; 441 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options")); 442 { 443 PetscBool flg; 444 PetscInt order; 445 CHKERRQ(TSBDFGetOrder(ts,&order)); 446 CHKERRQ(PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg)); 447 if (flg) CHKERRQ(TSBDFSetOrder(ts,order)); 448 } 449 CHKERRQ(PetscOptionsTail()); 450 PetscFunctionReturn(0); 451 } 452 453 static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer) 454 { 455 TS_BDF *bdf = (TS_BDF*)ts->data; 456 PetscBool iascii; 457 458 PetscFunctionBegin; 459 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 460 if (iascii) { 461 CHKERRQ(PetscViewerASCIIPrintf(viewer," Order=%D\n",bdf->order)); 462 } 463 PetscFunctionReturn(0); 464 } 465 466 /* ------------------------------------------------------------ */ 467 468 static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order) 469 { 470 TS_BDF *bdf = (TS_BDF*)ts->data; 471 472 PetscFunctionBegin; 473 if (order == bdf->order) PetscFunctionReturn(0); 474 PetscCheck(order >= 1 && order <= 6,PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order); 475 bdf->order = order; 476 PetscFunctionReturn(0); 477 } 478 479 static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order) 480 { 481 TS_BDF *bdf = (TS_BDF*)ts->data; 482 483 PetscFunctionBegin; 484 *order = bdf->order; 485 PetscFunctionReturn(0); 486 } 487 488 /* ------------------------------------------------------------ */ 489 490 /*MC 491 TSBDF - DAE solver using BDF methods 492 493 Level: beginner 494 495 .seealso: TS, TSCreate(), TSSetType() 496 M*/ 497 PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts) 498 { 499 TS_BDF *bdf; 500 501 PetscFunctionBegin; 502 ts->ops->reset = TSReset_BDF; 503 ts->ops->destroy = TSDestroy_BDF; 504 ts->ops->view = TSView_BDF; 505 ts->ops->setup = TSSetUp_BDF; 506 ts->ops->setfromoptions = TSSetFromOptions_BDF; 507 ts->ops->step = TSStep_BDF; 508 ts->ops->evaluatewlte = TSEvaluateWLTE_BDF; 509 ts->ops->rollback = TSRollBack_BDF; 510 ts->ops->interpolate = TSInterpolate_BDF; 511 ts->ops->snesfunction = SNESTSFormFunction_BDF; 512 ts->ops->snesjacobian = SNESTSFormJacobian_BDF; 513 ts->default_adapt_type = TSADAPTBASIC; 514 515 ts->usessnes = PETSC_TRUE; 516 517 CHKERRQ(PetscNewLog(ts,&bdf)); 518 ts->data = (void*)bdf; 519 520 bdf->status = TS_STEP_COMPLETE; 521 522 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF)); 523 CHKERRQ(PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF)); 524 CHKERRQ(TSBDFSetOrder(ts,2)); 525 PetscFunctionReturn(0); 526 } 527 528 /* ------------------------------------------------------------ */ 529 530 /*@ 531 TSBDFSetOrder - Set the order of the BDF method 532 533 Logically Collective on TS 534 535 Input Parameters: 536 + ts - timestepping context 537 - order - order of the method 538 539 Options Database: 540 . -ts_bdf_order <order> - select the order 541 542 Level: intermediate 543 544 @*/ 545 PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order) 546 { 547 PetscFunctionBegin; 548 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 549 PetscValidLogicalCollectiveInt(ts,order,2); 550 CHKERRQ(PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order))); 551 PetscFunctionReturn(0); 552 } 553 554 /*@ 555 TSBDFGetOrder - Get the order of the BDF method 556 557 Not Collective 558 559 Input Parameter: 560 . ts - timestepping context 561 562 Output Parameter: 563 . order - order of the method 564 565 Level: intermediate 566 567 @*/ 568 PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order) 569 { 570 PetscFunctionBegin; 571 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 572 PetscValidIntPointer(order,2); 573 CHKERRQ(PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order))); 574 PetscFunctionReturn(0); 575 } 576