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