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