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