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