1211a84d6SLisandro Dalcin /* 2211a84d6SLisandro Dalcin Code for timestepping with BDF methods 3211a84d6SLisandro Dalcin */ 4211a84d6SLisandro Dalcin #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 5211a84d6SLisandro Dalcin 6211a84d6SLisandro Dalcin static PetscBool cited = PETSC_FALSE; 7211a84d6SLisandro Dalcin static const char citation[] = 8211a84d6SLisandro Dalcin "@book{Brenan1995,\n" 9211a84d6SLisandro Dalcin " title = {Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations},\n" 10211a84d6SLisandro Dalcin " author = {Brenan, K. and Campbell, S. and Petzold, L.},\n" 11211a84d6SLisandro Dalcin " publisher = {Society for Industrial and Applied Mathematics},\n" 12211a84d6SLisandro Dalcin " year = {1995},\n" 13211a84d6SLisandro Dalcin " doi = {10.1137/1.9781611971224},\n}\n"; 14211a84d6SLisandro Dalcin 15211a84d6SLisandro Dalcin typedef struct { 16211a84d6SLisandro Dalcin PetscInt k,n; 17211a84d6SLisandro Dalcin PetscReal time[6+2]; 18211a84d6SLisandro Dalcin Vec work[6+2]; 19211a84d6SLisandro Dalcin PetscReal shift; 20211a84d6SLisandro Dalcin Vec vec_dot; 21211a84d6SLisandro Dalcin Vec vec_lte; 22211a84d6SLisandro Dalcin 23211a84d6SLisandro Dalcin PetscInt order; 24211a84d6SLisandro Dalcin PetscBool adapt; 25211a84d6SLisandro Dalcin TSStepStatus status; 26211a84d6SLisandro Dalcin } TS_BDF; 27211a84d6SLisandro Dalcin 28211a84d6SLisandro Dalcin 29211a84d6SLisandro Dalcin PETSC_STATIC_INLINE void LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[]) 30211a84d6SLisandro Dalcin { 31211a84d6SLisandro Dalcin PetscInt k,j; 32211a84d6SLisandro Dalcin for (k=0; k<n; k++) 33211a84d6SLisandro Dalcin for (L[k]=1, j=0; j<n; j++) 34211a84d6SLisandro Dalcin if (j != k) 35211a84d6SLisandro Dalcin L[k] *= (t - T[j])/(T[k] - T[j]); 36211a84d6SLisandro Dalcin } 37211a84d6SLisandro Dalcin 38211a84d6SLisandro Dalcin PETSC_STATIC_INLINE void LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[]) 39211a84d6SLisandro Dalcin { 40211a84d6SLisandro Dalcin PetscInt k,j,i; 41211a84d6SLisandro Dalcin for (k=0; k<n; k++) 42211a84d6SLisandro Dalcin for (dL[k]=0, j=0; j<n; j++) 43211a84d6SLisandro Dalcin if (j != k) { 44211a84d6SLisandro Dalcin PetscReal L = 1/(T[k] - T[j]); 45211a84d6SLisandro Dalcin for (i=0; i<n; i++) 46211a84d6SLisandro Dalcin if (i != j && i != k) 47211a84d6SLisandro Dalcin L *= (t - T[i])/(T[k] - T[i]); 48211a84d6SLisandro Dalcin dL[k] += L; 49211a84d6SLisandro Dalcin } 50211a84d6SLisandro Dalcin } 51211a84d6SLisandro Dalcin 52211a84d6SLisandro Dalcin #undef __FUNCT__ 53211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_Advance" 54211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_Advance(TS ts,PetscReal t,Vec X) 55211a84d6SLisandro Dalcin { 56211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 57211a84d6SLisandro Dalcin PetscInt i,n = (PetscInt)(sizeof(bdf->work)/sizeof(Vec)); 58211a84d6SLisandro Dalcin Vec tail = bdf->work[n-1]; 59211a84d6SLisandro Dalcin PetscErrorCode ierr; 60211a84d6SLisandro Dalcin 61211a84d6SLisandro Dalcin PetscFunctionBegin; 62211a84d6SLisandro Dalcin for (i=n-1; i>=2; i--) { 63211a84d6SLisandro Dalcin bdf->time[i] = bdf->time[i-1]; 64211a84d6SLisandro Dalcin bdf->work[i] = bdf->work[i-1]; 65211a84d6SLisandro Dalcin } 66211a84d6SLisandro Dalcin bdf->n = PetscMin(bdf->n+1,n-1); 67211a84d6SLisandro Dalcin bdf->time[1] = t; 68211a84d6SLisandro Dalcin bdf->work[1] = tail; 69211a84d6SLisandro Dalcin ierr = VecCopy(X,tail);CHKERRQ(ierr); 70211a84d6SLisandro Dalcin PetscFunctionReturn(0); 71211a84d6SLisandro Dalcin } 72211a84d6SLisandro Dalcin 73211a84d6SLisandro Dalcin #undef __FUNCT__ 74211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_VecDot" 75211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_VecDot(TS ts,PetscInt order,PetscReal t,Vec X,Vec Xdot,PetscReal *shift) 76211a84d6SLisandro Dalcin { 77211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 78211a84d6SLisandro Dalcin PetscInt i,n = order+1; 79211a84d6SLisandro Dalcin PetscReal time[7]; 80211a84d6SLisandro Dalcin Vec vecs[7]; 81211a84d6SLisandro Dalcin PetscScalar alpha[7]; 82211a84d6SLisandro Dalcin PetscErrorCode ierr; 83211a84d6SLisandro Dalcin 84211a84d6SLisandro Dalcin PetscFunctionBegin; 85211a84d6SLisandro Dalcin n = PetscMax(n,2); 86211a84d6SLisandro Dalcin time[0] = t; for (i=1; i<n; i++) time[i] = bdf->time[i]; 87211a84d6SLisandro Dalcin vecs[0] = X; for (i=1; i<n; i++) vecs[i] = bdf->work[i]; 88211a84d6SLisandro Dalcin LagrangeBasisDers(n,t,time,alpha); 89211a84d6SLisandro Dalcin ierr = VecZeroEntries(Xdot);CHKERRQ(ierr); 90211a84d6SLisandro Dalcin ierr = VecMAXPY(Xdot,n,alpha,vecs);CHKERRQ(ierr); 91211a84d6SLisandro Dalcin if (shift) *shift = PetscRealPart(alpha[0]); 92211a84d6SLisandro Dalcin PetscFunctionReturn(0); 93211a84d6SLisandro Dalcin } 94211a84d6SLisandro Dalcin 95211a84d6SLisandro Dalcin #undef __FUNCT__ 96211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_VecLTE" 97211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_VecLTE(TS ts,PetscInt order,Vec lte) 98211a84d6SLisandro Dalcin { 99211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 100211a84d6SLisandro Dalcin PetscInt i,n = order+1; 101211a84d6SLisandro Dalcin PetscReal *time = bdf->time; 102211a84d6SLisandro Dalcin Vec *vecs = bdf->work; 103211a84d6SLisandro Dalcin PetscScalar a[8],b[8],alpha[8]; 104211a84d6SLisandro Dalcin PetscErrorCode ierr; 105211a84d6SLisandro Dalcin 106211a84d6SLisandro Dalcin PetscFunctionBegin; 107211a84d6SLisandro Dalcin LagrangeBasisDers(n+0,time[0],time,a); a[n] =0; 108211a84d6SLisandro Dalcin LagrangeBasisDers(n+1,time[0],time,b); 109211a84d6SLisandro Dalcin for (i=0; i<n+1; i++) alpha[i] = (a[i]-b[i])/a[0]; 110211a84d6SLisandro Dalcin ierr = VecZeroEntries(lte);CHKERRQ(ierr); 111211a84d6SLisandro Dalcin ierr = VecMAXPY(lte,n+1,alpha,vecs);CHKERRQ(ierr); 112211a84d6SLisandro Dalcin PetscFunctionReturn(0); 113211a84d6SLisandro Dalcin } 114211a84d6SLisandro Dalcin 115211a84d6SLisandro Dalcin #undef __FUNCT__ 116211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_Extrapolate" 117211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_Extrapolate(TS ts,PetscInt order,PetscReal t,Vec X) 118211a84d6SLisandro Dalcin { 119211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 120211a84d6SLisandro Dalcin PetscInt n = order+1; 121211a84d6SLisandro Dalcin PetscReal *time = bdf->time+1; 122211a84d6SLisandro Dalcin Vec *vecs = bdf->work+1; 123211a84d6SLisandro Dalcin PetscScalar alpha[7]; 124211a84d6SLisandro Dalcin PetscErrorCode ierr; 125211a84d6SLisandro Dalcin 126211a84d6SLisandro Dalcin PetscFunctionBegin; 127211a84d6SLisandro Dalcin n = PetscMin(n,bdf->n); 128211a84d6SLisandro Dalcin LagrangeBasisVals(n,t,time,alpha); 129211a84d6SLisandro Dalcin ierr = VecZeroEntries(X);CHKERRQ(ierr); 130211a84d6SLisandro Dalcin ierr = VecMAXPY(X,n,alpha,vecs);CHKERRQ(ierr); 131211a84d6SLisandro Dalcin PetscFunctionReturn(0); 132211a84d6SLisandro Dalcin } 133211a84d6SLisandro Dalcin 134211a84d6SLisandro Dalcin #undef __FUNCT__ 135211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_Interpolate" 136211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_Interpolate(TS ts,PetscInt order,PetscReal t,Vec X) 137211a84d6SLisandro Dalcin { 138211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 139211a84d6SLisandro Dalcin PetscInt n = order+1; 140211a84d6SLisandro Dalcin PetscReal *time = bdf->time; 141211a84d6SLisandro Dalcin Vec *vecs = bdf->work; 142211a84d6SLisandro Dalcin PetscScalar alpha[7]; 143211a84d6SLisandro Dalcin PetscErrorCode ierr; 144211a84d6SLisandro Dalcin 145211a84d6SLisandro Dalcin PetscFunctionBegin; 146211a84d6SLisandro Dalcin LagrangeBasisVals(n,t,time,alpha); 147211a84d6SLisandro Dalcin ierr = VecZeroEntries(X);CHKERRQ(ierr); 148211a84d6SLisandro Dalcin ierr = VecMAXPY(X,n,alpha,vecs);CHKERRQ(ierr); 149211a84d6SLisandro Dalcin PetscFunctionReturn(0); 150211a84d6SLisandro Dalcin } 151211a84d6SLisandro Dalcin 152211a84d6SLisandro Dalcin #undef __FUNCT__ 153211a84d6SLisandro Dalcin #define __FUNCT__ "TS_SNESSolve" 154211a84d6SLisandro Dalcin static PetscErrorCode TS_SNESSolve(TS ts,Vec b,Vec x) 155211a84d6SLisandro Dalcin { 156211a84d6SLisandro Dalcin PetscInt nits,lits; 157211a84d6SLisandro Dalcin PetscErrorCode ierr; 158211a84d6SLisandro Dalcin 159211a84d6SLisandro Dalcin PetscFunctionBegin; 160211a84d6SLisandro Dalcin ierr = SNESSolve(ts->snes,b,x);CHKERRQ(ierr); 161211a84d6SLisandro Dalcin ierr = SNESGetIterationNumber(ts->snes,&nits);CHKERRQ(ierr); 162211a84d6SLisandro Dalcin ierr = SNESGetLinearSolveIterations(ts->snes,&lits);CHKERRQ(ierr); 163211a84d6SLisandro Dalcin ts->snes_its += nits; ts->ksp_its += lits; 164211a84d6SLisandro Dalcin PetscFunctionReturn(0); 165211a84d6SLisandro Dalcin } 166211a84d6SLisandro Dalcin 167211a84d6SLisandro Dalcin #undef __FUNCT__ 168211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDF_Restart" 169211a84d6SLisandro Dalcin static PetscErrorCode TSBDF_Restart(TS ts,PetscBool *accept) 170211a84d6SLisandro Dalcin { 171211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 172211a84d6SLisandro Dalcin PetscErrorCode ierr; 173211a84d6SLisandro Dalcin 174211a84d6SLisandro Dalcin PetscFunctionBegin; 175211a84d6SLisandro Dalcin bdf->k = 1; bdf->n = 0; 176211a84d6SLisandro Dalcin ierr = TSBDF_Advance(ts,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 177211a84d6SLisandro Dalcin 178211a84d6SLisandro Dalcin bdf->time[0] = ts->ptime + ts->time_step/2; 179211a84d6SLisandro Dalcin ierr = VecCopy(bdf->work[1],bdf->work[0]);CHKERRQ(ierr); 180211a84d6SLisandro Dalcin ierr = TSPreStage(ts,bdf->time[0]);CHKERRQ(ierr); 181211a84d6SLisandro Dalcin ierr = TS_SNESSolve(ts,NULL,bdf->work[0]);CHKERRQ(ierr); 182211a84d6SLisandro Dalcin ierr = TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);CHKERRQ(ierr); 183211a84d6SLisandro Dalcin ierr = TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],accept);CHKERRQ(ierr); 184211a84d6SLisandro Dalcin if (!*accept) PetscFunctionReturn(0); 185211a84d6SLisandro Dalcin 186211a84d6SLisandro Dalcin bdf->k = PetscMin(2,bdf->order); bdf->n++; 187211a84d6SLisandro Dalcin ierr = VecCopy(bdf->work[0],bdf->work[2]);CHKERRQ(ierr); 188211a84d6SLisandro Dalcin bdf->time[2] = bdf->time[0]; 189211a84d6SLisandro Dalcin PetscFunctionReturn(0); 190211a84d6SLisandro Dalcin } 191211a84d6SLisandro Dalcin 192211a84d6SLisandro Dalcin static const char *const BDF_SchemeName[] = {"", "1", "2", "3", "4", "5", "6"}; 193211a84d6SLisandro Dalcin 194211a84d6SLisandro Dalcin #undef __FUNCT__ 195211a84d6SLisandro Dalcin #define __FUNCT__ "TSStep_BDF" 196211a84d6SLisandro Dalcin static PetscErrorCode TSStep_BDF(TS ts) 197211a84d6SLisandro Dalcin { 198211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 199211a84d6SLisandro Dalcin PetscInt rejections = 0; 200211a84d6SLisandro Dalcin PetscBool stageok,accept = PETSC_TRUE; 201211a84d6SLisandro Dalcin PetscReal next_time_step = ts->time_step; 202211a84d6SLisandro Dalcin PetscErrorCode ierr; 203211a84d6SLisandro Dalcin 204211a84d6SLisandro Dalcin PetscFunctionBegin; 205211a84d6SLisandro Dalcin ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 206211a84d6SLisandro Dalcin 207211a84d6SLisandro Dalcin if (!ts->steprollback && !ts->steprestart) { 208211a84d6SLisandro Dalcin bdf->k = PetscMin(bdf->k+1,bdf->order); 209211a84d6SLisandro Dalcin ierr = TSBDF_Advance(ts,ts->ptime,ts->vec_sol);CHKERRQ(ierr); 210211a84d6SLisandro Dalcin } 211211a84d6SLisandro Dalcin 212211a84d6SLisandro Dalcin bdf->status = TS_STEP_INCOMPLETE; 213211a84d6SLisandro Dalcin while (!ts->reason && bdf->status != TS_STEP_COMPLETE) { 214211a84d6SLisandro Dalcin 215211a84d6SLisandro Dalcin if (ts->steprestart) { 216211a84d6SLisandro Dalcin ierr = TSBDF_Restart(ts,&stageok);CHKERRQ(ierr); 217211a84d6SLisandro Dalcin if (!stageok) goto reject_step; 218211a84d6SLisandro Dalcin } 219211a84d6SLisandro Dalcin 220211a84d6SLisandro Dalcin bdf->time[0] = ts->ptime + ts->time_step; 221211a84d6SLisandro Dalcin ierr = TSBDF_Extrapolate(ts,bdf->k-(accept?0:1),bdf->time[0],bdf->work[0]);CHKERRQ(ierr); 222211a84d6SLisandro Dalcin ierr = TSPreStage(ts,bdf->time[0]);CHKERRQ(ierr); 223211a84d6SLisandro Dalcin ierr = TS_SNESSolve(ts,NULL,bdf->work[0]);CHKERRQ(ierr); 224211a84d6SLisandro Dalcin ierr = TSPostStage(ts,bdf->time[0],0,&bdf->work[0]);CHKERRQ(ierr); 225211a84d6SLisandro Dalcin ierr = TSAdaptCheckStage(ts->adapt,ts,bdf->time[0],bdf->work[0],&stageok);CHKERRQ(ierr); 226211a84d6SLisandro Dalcin if (!stageok) goto reject_step; 227211a84d6SLisandro Dalcin 228211a84d6SLisandro Dalcin bdf->status = TS_STEP_PENDING; 229211a84d6SLisandro Dalcin ierr = TSAdaptCandidatesClear(ts->adapt);CHKERRQ(ierr); 230211a84d6SLisandro Dalcin ierr = TSAdaptCandidateAdd(ts->adapt,BDF_SchemeName[bdf->k],bdf->k,1,1.0,1.0,PETSC_TRUE);CHKERRQ(ierr); 231211a84d6SLisandro Dalcin ierr = TSAdaptChoose(ts->adapt,ts,ts->time_step,NULL,&next_time_step,&accept);CHKERRQ(ierr); 232211a84d6SLisandro Dalcin bdf->status = accept ? TS_STEP_COMPLETE : TS_STEP_INCOMPLETE; 233211a84d6SLisandro Dalcin if (!accept) { ts->time_step = next_time_step; goto reject_step; } 234211a84d6SLisandro Dalcin 235211a84d6SLisandro Dalcin ierr = VecCopy(bdf->work[0],ts->vec_sol);CHKERRQ(ierr); 236211a84d6SLisandro Dalcin ts->ptime += ts->time_step; 237211a84d6SLisandro Dalcin ts->time_step = next_time_step; 238211a84d6SLisandro Dalcin break; 239211a84d6SLisandro Dalcin 240211a84d6SLisandro Dalcin reject_step: 241211a84d6SLisandro Dalcin ts->reject++; accept = PETSC_FALSE; 242211a84d6SLisandro Dalcin if (!ts->reason && ++rejections > ts->max_reject && ts->max_reject >= 0) { 243211a84d6SLisandro Dalcin ierr = PetscInfo2(ts,"Step=%D, step rejections %D greater than current TS allowed, stopping solve\n",ts->steps,rejections);CHKERRQ(ierr); 244211a84d6SLisandro Dalcin ts->reason = TS_DIVERGED_STEP_REJECTED; 245211a84d6SLisandro Dalcin } 246211a84d6SLisandro Dalcin } 247211a84d6SLisandro Dalcin PetscFunctionReturn(0); 248211a84d6SLisandro Dalcin } 249211a84d6SLisandro Dalcin 250211a84d6SLisandro Dalcin #undef __FUNCT__ 251211a84d6SLisandro Dalcin #define __FUNCT__ "TSInterpolate_BDF" 252211a84d6SLisandro Dalcin static PetscErrorCode TSInterpolate_BDF(TS ts,PetscReal t,Vec X) 253211a84d6SLisandro Dalcin { 254211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 255211a84d6SLisandro Dalcin PetscErrorCode ierr; 256211a84d6SLisandro Dalcin 257211a84d6SLisandro Dalcin PetscFunctionBegin; 258211a84d6SLisandro Dalcin ierr = TSBDF_Interpolate(ts,bdf->k,t,X);CHKERRQ(ierr); 259211a84d6SLisandro Dalcin PetscFunctionReturn(0); 260211a84d6SLisandro Dalcin } 261211a84d6SLisandro Dalcin 262211a84d6SLisandro Dalcin #undef __FUNCT__ 263211a84d6SLisandro Dalcin #define __FUNCT__ "TSEvaluateWLTE_BDF" 264211a84d6SLisandro Dalcin static PetscErrorCode TSEvaluateWLTE_BDF(TS ts,NormType wnormtype,PetscInt *order,PetscReal *wlte) 265211a84d6SLisandro Dalcin { 266211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 267211a84d6SLisandro Dalcin PetscInt k = bdf->k; 268*7453f775SEmil Constantinescu PetscReal wltea,wlter; 269211a84d6SLisandro Dalcin Vec X = bdf->work[0], Y = bdf->vec_lte; 270211a84d6SLisandro Dalcin PetscErrorCode ierr; 271211a84d6SLisandro Dalcin 272211a84d6SLisandro Dalcin PetscFunctionBegin; 273211a84d6SLisandro Dalcin k = PetscMin(k,bdf->n-1); 274211a84d6SLisandro Dalcin ierr = TSBDF_VecLTE(ts,k,Y);CHKERRQ(ierr); 275211a84d6SLisandro Dalcin ierr = VecAXPY(Y,1,X);CHKERRQ(ierr); 276*7453f775SEmil Constantinescu ierr = TSErrorWeightedNorm(ts,X,Y,wnormtype,wlte,&wltea,&wlter);CHKERRQ(ierr); 277211a84d6SLisandro Dalcin if (order) *order = k + 1; 278211a84d6SLisandro Dalcin PetscFunctionReturn(0); 279211a84d6SLisandro Dalcin } 280211a84d6SLisandro Dalcin 281211a84d6SLisandro Dalcin #undef __FUNCT__ 282211a84d6SLisandro Dalcin #define __FUNCT__ "TSRollBack_BDF" 283211a84d6SLisandro Dalcin static PetscErrorCode TSRollBack_BDF(TS ts) 284211a84d6SLisandro Dalcin { 285211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 286211a84d6SLisandro Dalcin PetscErrorCode ierr; 287211a84d6SLisandro Dalcin 288211a84d6SLisandro Dalcin PetscFunctionBegin; 289211a84d6SLisandro Dalcin ierr = VecCopy(bdf->work[1],ts->vec_sol);CHKERRQ(ierr); 290211a84d6SLisandro Dalcin PetscFunctionReturn(0); 291211a84d6SLisandro Dalcin } 292211a84d6SLisandro Dalcin 293211a84d6SLisandro Dalcin #undef __FUNCT__ 294211a84d6SLisandro Dalcin #define __FUNCT__ "SNESTSFormFunction_BDF" 295211a84d6SLisandro Dalcin static PetscErrorCode SNESTSFormFunction_BDF(PETSC_UNUSED SNES snes,Vec X,Vec F,TS ts) 296211a84d6SLisandro Dalcin { 297211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 298211a84d6SLisandro Dalcin PetscInt k = bdf->k; 299211a84d6SLisandro Dalcin PetscReal t = bdf->time[0]; 300211a84d6SLisandro Dalcin Vec V = bdf->vec_dot; 301211a84d6SLisandro Dalcin PetscErrorCode ierr; 302211a84d6SLisandro Dalcin 303211a84d6SLisandro Dalcin PetscFunctionBegin; 304211a84d6SLisandro Dalcin ierr = TSBDF_VecDot(ts,k,t,X,V,&bdf->shift);CHKERRQ(ierr); 305211a84d6SLisandro Dalcin /* F = Function(t,X,V) */ 306211a84d6SLisandro Dalcin ierr = TSComputeIFunction(ts,t,X,V,F,PETSC_FALSE);CHKERRQ(ierr); 307211a84d6SLisandro Dalcin PetscFunctionReturn(0); 308211a84d6SLisandro Dalcin } 309211a84d6SLisandro Dalcin 310211a84d6SLisandro Dalcin #undef __FUNCT__ 311211a84d6SLisandro Dalcin #define __FUNCT__ "SNESTSFormJacobian_BDF" 312211a84d6SLisandro Dalcin static PetscErrorCode SNESTSFormJacobian_BDF(PETSC_UNUSED SNES snes, 313211a84d6SLisandro Dalcin PETSC_UNUSED Vec X, 314211a84d6SLisandro Dalcin Mat J,Mat P, 315211a84d6SLisandro Dalcin TS ts) 316211a84d6SLisandro Dalcin { 317211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 318211a84d6SLisandro Dalcin PetscReal t = bdf->time[0]; 319211a84d6SLisandro Dalcin Vec V = bdf->vec_dot; 320211a84d6SLisandro Dalcin PetscReal dVdX = bdf->shift; 321211a84d6SLisandro Dalcin PetscErrorCode ierr; 322211a84d6SLisandro Dalcin 323211a84d6SLisandro Dalcin PetscFunctionBegin; 324211a84d6SLisandro Dalcin /* J,P = Jacobian(t,X,V) */ 325211a84d6SLisandro Dalcin ierr = TSComputeIJacobian(ts,t,X,V,dVdX,J,P,PETSC_FALSE);CHKERRQ(ierr); 326211a84d6SLisandro Dalcin PetscFunctionReturn(0); 327211a84d6SLisandro Dalcin } 328211a84d6SLisandro Dalcin 329211a84d6SLisandro Dalcin #undef __FUNCT__ 330211a84d6SLisandro Dalcin #define __FUNCT__ "TSReset_BDF" 331211a84d6SLisandro Dalcin static PetscErrorCode TSReset_BDF(TS ts) 332211a84d6SLisandro Dalcin { 333211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 334211a84d6SLisandro Dalcin size_t i,n = sizeof(bdf->work)/sizeof(Vec); 335211a84d6SLisandro Dalcin PetscErrorCode ierr; 336211a84d6SLisandro Dalcin 337211a84d6SLisandro Dalcin PetscFunctionBegin; 338211a84d6SLisandro Dalcin for (i=0; i<n; i++) {ierr = VecDestroy(&bdf->work[i]);CHKERRQ(ierr);} 339211a84d6SLisandro Dalcin ierr = VecDestroy(&bdf->vec_dot);CHKERRQ(ierr); 340211a84d6SLisandro Dalcin ierr = VecDestroy(&bdf->vec_lte);CHKERRQ(ierr); 341211a84d6SLisandro Dalcin PetscFunctionReturn(0); 342211a84d6SLisandro Dalcin } 343211a84d6SLisandro Dalcin 344211a84d6SLisandro Dalcin #undef __FUNCT__ 345211a84d6SLisandro Dalcin #define __FUNCT__ "TSDestroy_BDF" 346211a84d6SLisandro Dalcin static PetscErrorCode TSDestroy_BDF(TS ts) 347211a84d6SLisandro Dalcin { 348211a84d6SLisandro Dalcin PetscErrorCode ierr; 349211a84d6SLisandro Dalcin 350211a84d6SLisandro Dalcin PetscFunctionBegin; 351211a84d6SLisandro Dalcin ierr = TSReset_BDF(ts);CHKERRQ(ierr); 352211a84d6SLisandro Dalcin ierr = PetscFree(ts->data);CHKERRQ(ierr); 353211a84d6SLisandro Dalcin ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",NULL);CHKERRQ(ierr); 354211a84d6SLisandro Dalcin ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",NULL);CHKERRQ(ierr); 355211a84d6SLisandro Dalcin ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFUseAdapt_C",NULL);CHKERRQ(ierr); 356211a84d6SLisandro Dalcin PetscFunctionReturn(0); 357211a84d6SLisandro Dalcin } 358211a84d6SLisandro Dalcin 359211a84d6SLisandro Dalcin #undef __FUNCT__ 360211a84d6SLisandro Dalcin #define __FUNCT__ "TSSetUp_BDF" 361211a84d6SLisandro Dalcin static PetscErrorCode TSSetUp_BDF(TS ts) 362211a84d6SLisandro Dalcin { 363211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 364211a84d6SLisandro Dalcin size_t i,n = sizeof(bdf->work)/sizeof(Vec); 365211a84d6SLisandro Dalcin PetscErrorCode ierr; 366211a84d6SLisandro Dalcin 367211a84d6SLisandro Dalcin PetscFunctionBegin; 368211a84d6SLisandro Dalcin bdf->k = bdf->n = 0; 369211a84d6SLisandro Dalcin for (i=0; i<n; i++) {ierr = VecDuplicate(ts->vec_sol,&bdf->work[i]);CHKERRQ(ierr);} 370211a84d6SLisandro Dalcin ierr = VecDuplicate(ts->vec_sol,&bdf->vec_dot);CHKERRQ(ierr); 371211a84d6SLisandro Dalcin ierr = VecDuplicate(ts->vec_sol,&bdf->vec_lte);CHKERRQ(ierr); 372211a84d6SLisandro Dalcin 373211a84d6SLisandro Dalcin ierr = TSGetAdapt(ts,&ts->adapt);CHKERRQ(ierr); 374211a84d6SLisandro Dalcin ierr = TSAdaptCandidatesClear(ts->adapt);CHKERRQ(ierr); 375211a84d6SLisandro Dalcin if (!bdf->adapt) { 376211a84d6SLisandro Dalcin ierr = TSAdaptSetType(ts->adapt,TSADAPTNONE);CHKERRQ(ierr); 377211a84d6SLisandro Dalcin } else { 378211a84d6SLisandro Dalcin PetscReal low,high; 379211a84d6SLisandro Dalcin ierr = TSAdaptBasicGetClip(ts->adapt,&low,&high);CHKERRQ(ierr); 380211a84d6SLisandro Dalcin high = PetscMin(high,2.0); 381211a84d6SLisandro Dalcin ierr = TSAdaptBasicSetClip(ts->adapt,low,high);CHKERRQ(ierr); 382211a84d6SLisandro Dalcin if (ts->exact_final_time == TS_EXACTFINALTIME_UNSPECIFIED) 383211a84d6SLisandro Dalcin ts->exact_final_time = TS_EXACTFINALTIME_MATCHSTEP; 384211a84d6SLisandro Dalcin } 385211a84d6SLisandro Dalcin 386211a84d6SLisandro Dalcin ierr = TSGetSNES(ts,&ts->snes);CHKERRQ(ierr); 387211a84d6SLisandro Dalcin PetscFunctionReturn(0); 388211a84d6SLisandro Dalcin } 389211a84d6SLisandro Dalcin 390211a84d6SLisandro Dalcin #undef __FUNCT__ 391211a84d6SLisandro Dalcin #define __FUNCT__ "TSSetFromOptions_BDF" 392211a84d6SLisandro Dalcin static PetscErrorCode TSSetFromOptions_BDF(PetscOptionItems *PetscOptionsObject,TS ts) 393211a84d6SLisandro Dalcin { 394211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 395211a84d6SLisandro Dalcin PetscErrorCode ierr; 396211a84d6SLisandro Dalcin 397211a84d6SLisandro Dalcin PetscFunctionBegin; 398211a84d6SLisandro Dalcin ierr = PetscOptionsHead(PetscOptionsObject,"BDF ODE solver options");CHKERRQ(ierr); 399211a84d6SLisandro Dalcin { 400211a84d6SLisandro Dalcin PetscBool flg; 401211a84d6SLisandro Dalcin PetscInt order = bdf->order; 402211a84d6SLisandro Dalcin PetscBool adapt = bdf->adapt; 403211a84d6SLisandro Dalcin ierr = PetscOptionsInt("-ts_bdf_order","Order of the BDF method","TSBDFSetOrder",order,&order,&flg);CHKERRQ(ierr); 404211a84d6SLisandro Dalcin if (flg) {ierr = TSBDFSetOrder(ts,order);CHKERRQ(ierr);} 405211a84d6SLisandro Dalcin ierr = PetscOptionsBool("-ts_bdf_adapt","Use time-step adaptivity with the BDF method","TSBDFUseAdapt",adapt,&adapt,&flg);CHKERRQ(ierr); 406211a84d6SLisandro Dalcin if (flg) {ierr = TSBDFUseAdapt(ts,adapt);CHKERRQ(ierr);} 407211a84d6SLisandro Dalcin } 408211a84d6SLisandro Dalcin ierr = PetscOptionsTail();CHKERRQ(ierr); 409211a84d6SLisandro Dalcin PetscFunctionReturn(0); 410211a84d6SLisandro Dalcin } 411211a84d6SLisandro Dalcin 412211a84d6SLisandro Dalcin #undef __FUNCT__ 413211a84d6SLisandro Dalcin #define __FUNCT__ "TSView_BDF" 414211a84d6SLisandro Dalcin static PetscErrorCode TSView_BDF(TS ts,PetscViewer viewer) 415211a84d6SLisandro Dalcin { 416211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 417211a84d6SLisandro Dalcin PetscBool iascii; 418211a84d6SLisandro Dalcin PetscErrorCode ierr; 419211a84d6SLisandro Dalcin 420211a84d6SLisandro Dalcin PetscFunctionBegin; 421211a84d6SLisandro Dalcin ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 422211a84d6SLisandro Dalcin if (iascii) { 423211a84d6SLisandro Dalcin ierr = PetscViewerASCIIPrintf(viewer," Order=%D\n",bdf->order);CHKERRQ(ierr); 424211a84d6SLisandro Dalcin } 425211a84d6SLisandro Dalcin if (ts->adapt) {ierr = TSAdaptView(ts->adapt,viewer);CHKERRQ(ierr);} 426211a84d6SLisandro Dalcin if (ts->snes) {ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr);} 427211a84d6SLisandro Dalcin PetscFunctionReturn(0); 428211a84d6SLisandro Dalcin } 429211a84d6SLisandro Dalcin 430211a84d6SLisandro Dalcin /* ------------------------------------------------------------ */ 431211a84d6SLisandro Dalcin 432211a84d6SLisandro Dalcin #undef __FUNCT__ 433211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFSetOrder_BDF" 434211a84d6SLisandro Dalcin static PetscErrorCode TSBDFSetOrder_BDF(TS ts,PetscInt order) 435211a84d6SLisandro Dalcin { 436211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 437211a84d6SLisandro Dalcin 438211a84d6SLisandro Dalcin PetscFunctionBegin; 439211a84d6SLisandro Dalcin if (order == bdf->order) PetscFunctionReturn(0); 440211a84d6SLisandro Dalcin if (order < 1 || order > 6) SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_OUTOFRANGE,"BDF Order %D not implemented",order); 441211a84d6SLisandro Dalcin bdf->order = order; 442211a84d6SLisandro Dalcin PetscFunctionReturn(0); 443211a84d6SLisandro Dalcin } 444211a84d6SLisandro Dalcin 445211a84d6SLisandro Dalcin #undef __FUNCT__ 446211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFGetOrder_BDF" 447211a84d6SLisandro Dalcin static PetscErrorCode TSBDFGetOrder_BDF(TS ts,PetscInt *order) 448211a84d6SLisandro Dalcin { 449211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 450211a84d6SLisandro Dalcin 451211a84d6SLisandro Dalcin PetscFunctionBegin; 452211a84d6SLisandro Dalcin *order = bdf->order; 453211a84d6SLisandro Dalcin PetscFunctionReturn(0); 454211a84d6SLisandro Dalcin } 455211a84d6SLisandro Dalcin 456211a84d6SLisandro Dalcin #undef __FUNCT__ 457211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFUseAdapt_BDF" 458211a84d6SLisandro Dalcin static PetscErrorCode TSBDFUseAdapt_BDF(TS ts,PetscBool use) 459211a84d6SLisandro Dalcin { 460211a84d6SLisandro Dalcin TS_BDF *bdf = (TS_BDF*)ts->data; 461211a84d6SLisandro Dalcin 462211a84d6SLisandro Dalcin PetscFunctionBegin; 463211a84d6SLisandro Dalcin if (use == bdf->adapt) PetscFunctionReturn(0); 464211a84d6SLisandro Dalcin if (ts->setupcalled) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_ORDER,"Cannot change adaptivity after TSSetUp()"); 465211a84d6SLisandro Dalcin bdf->adapt = use; 466211a84d6SLisandro Dalcin PetscFunctionReturn(0); 467211a84d6SLisandro Dalcin } 468211a84d6SLisandro Dalcin 469211a84d6SLisandro Dalcin /* ------------------------------------------------------------ */ 470211a84d6SLisandro Dalcin 471211a84d6SLisandro Dalcin /*MC 472211a84d6SLisandro Dalcin TSBDF - DAE solver using BDF methods 473211a84d6SLisandro Dalcin 474211a84d6SLisandro Dalcin Level: beginner 475211a84d6SLisandro Dalcin 476211a84d6SLisandro Dalcin .seealso: TS, TSCreate(), TSSetType() 477211a84d6SLisandro Dalcin M*/ 478211a84d6SLisandro Dalcin #undef __FUNCT__ 479211a84d6SLisandro Dalcin #define __FUNCT__ "TSCreate_BDF" 480211a84d6SLisandro Dalcin PETSC_EXTERN PetscErrorCode TSCreate_BDF(TS ts) 481211a84d6SLisandro Dalcin { 482211a84d6SLisandro Dalcin TS_BDF *bdf; 483211a84d6SLisandro Dalcin PetscErrorCode ierr; 484211a84d6SLisandro Dalcin 485211a84d6SLisandro Dalcin PetscFunctionBegin; 486211a84d6SLisandro Dalcin ts->ops->reset = TSReset_BDF; 487211a84d6SLisandro Dalcin ts->ops->destroy = TSDestroy_BDF; 488211a84d6SLisandro Dalcin ts->ops->view = TSView_BDF; 489211a84d6SLisandro Dalcin ts->ops->setup = TSSetUp_BDF; 490211a84d6SLisandro Dalcin ts->ops->setfromoptions = TSSetFromOptions_BDF; 491211a84d6SLisandro Dalcin ts->ops->step = TSStep_BDF; 492211a84d6SLisandro Dalcin ts->ops->evaluatewlte = TSEvaluateWLTE_BDF; 493211a84d6SLisandro Dalcin ts->ops->rollback = TSRollBack_BDF; 494211a84d6SLisandro Dalcin ts->ops->interpolate = TSInterpolate_BDF; 495211a84d6SLisandro Dalcin ts->ops->snesfunction = SNESTSFormFunction_BDF; 496211a84d6SLisandro Dalcin ts->ops->snesjacobian = SNESTSFormJacobian_BDF; 497211a84d6SLisandro Dalcin 498211a84d6SLisandro Dalcin ierr = PetscNewLog(ts,&bdf);CHKERRQ(ierr); 499211a84d6SLisandro Dalcin ts->data = (void*)bdf; 500211a84d6SLisandro Dalcin 501211a84d6SLisandro Dalcin bdf->order = 2; 502211a84d6SLisandro Dalcin bdf->adapt = PETSC_FALSE; 503211a84d6SLisandro Dalcin bdf->status = TS_STEP_COMPLETE; 504211a84d6SLisandro Dalcin 505211a84d6SLisandro Dalcin ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFSetOrder_C",TSBDFSetOrder_BDF);CHKERRQ(ierr); 506211a84d6SLisandro Dalcin ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFGetOrder_C",TSBDFGetOrder_BDF);CHKERRQ(ierr); 507211a84d6SLisandro Dalcin ierr = PetscObjectComposeFunction((PetscObject)ts,"TSBDFUseAdapt_C",TSBDFUseAdapt_BDF);CHKERRQ(ierr); 508211a84d6SLisandro Dalcin PetscFunctionReturn(0); 509211a84d6SLisandro Dalcin } 510211a84d6SLisandro Dalcin 511211a84d6SLisandro Dalcin /* ------------------------------------------------------------ */ 512211a84d6SLisandro Dalcin 513211a84d6SLisandro Dalcin #undef __FUNCT__ 514211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFSetOrder" 515211a84d6SLisandro Dalcin /*@ 516211a84d6SLisandro Dalcin TSBDFSetOrder - Set the order of the BDF method 517211a84d6SLisandro Dalcin 518211a84d6SLisandro Dalcin Logically Collective on TS 519211a84d6SLisandro Dalcin 520211a84d6SLisandro Dalcin Input Parameter: 521211a84d6SLisandro Dalcin + ts - timestepping context 522211a84d6SLisandro Dalcin - order - order of the method 523211a84d6SLisandro Dalcin 524211a84d6SLisandro Dalcin Options Database: 525211a84d6SLisandro Dalcin . -ts_bdf_order <order> 526211a84d6SLisandro Dalcin 527211a84d6SLisandro Dalcin Level: intermediate 528211a84d6SLisandro Dalcin 529211a84d6SLisandro Dalcin @*/ 530211a84d6SLisandro Dalcin PetscErrorCode TSBDFSetOrder(TS ts,PetscInt order) 531211a84d6SLisandro Dalcin { 532211a84d6SLisandro Dalcin PetscErrorCode ierr; 533211a84d6SLisandro Dalcin 534211a84d6SLisandro Dalcin PetscFunctionBegin; 535211a84d6SLisandro Dalcin PetscValidHeaderSpecific(ts,TS_CLASSID,1); 536211a84d6SLisandro Dalcin PetscValidLogicalCollectiveInt(ts,order,2); 537211a84d6SLisandro Dalcin ierr = PetscTryMethod(ts,"TSBDFSetOrder_C",(TS,PetscInt),(ts,order));CHKERRQ(ierr); 538211a84d6SLisandro Dalcin PetscFunctionReturn(0); 539211a84d6SLisandro Dalcin } 540211a84d6SLisandro Dalcin 541211a84d6SLisandro Dalcin #undef __FUNCT__ 542211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFGetOrder" 543211a84d6SLisandro Dalcin /*@ 544211a84d6SLisandro Dalcin TSBDFGetOrder - Get the order of the BDF method 545211a84d6SLisandro Dalcin 546211a84d6SLisandro Dalcin Not Collective 547211a84d6SLisandro Dalcin 548211a84d6SLisandro Dalcin Input Parameter: 549211a84d6SLisandro Dalcin . ts - timestepping context 550211a84d6SLisandro Dalcin 551211a84d6SLisandro Dalcin Output Parameter: 552211a84d6SLisandro Dalcin . order - order of the method 553211a84d6SLisandro Dalcin 554211a84d6SLisandro Dalcin Level: intermediate 555211a84d6SLisandro Dalcin 556211a84d6SLisandro Dalcin @*/ 557211a84d6SLisandro Dalcin PetscErrorCode TSBDFGetOrder(TS ts,PetscInt *order) 558211a84d6SLisandro Dalcin { 559211a84d6SLisandro Dalcin PetscErrorCode ierr; 560211a84d6SLisandro Dalcin 561211a84d6SLisandro Dalcin PetscFunctionBegin; 562211a84d6SLisandro Dalcin PetscValidHeaderSpecific(ts,TS_CLASSID,1); 563211a84d6SLisandro Dalcin PetscValidIntPointer(order,2); 564211a84d6SLisandro Dalcin ierr = PetscUseMethod(ts,"TSBDFGetOrder_C",(TS,PetscInt*),(ts,order));CHKERRQ(ierr); 565211a84d6SLisandro Dalcin PetscFunctionReturn(0); 566211a84d6SLisandro Dalcin } 567211a84d6SLisandro Dalcin 568211a84d6SLisandro Dalcin #undef __FUNCT__ 569211a84d6SLisandro Dalcin #define __FUNCT__ "TSBDFUseAdapt" 570211a84d6SLisandro Dalcin /*@ 571211a84d6SLisandro Dalcin TSBDFUseAdapt - Use time-step adaptivity with the BDF method 572211a84d6SLisandro Dalcin 573211a84d6SLisandro Dalcin Logically Collective on TS 574211a84d6SLisandro Dalcin 575211a84d6SLisandro Dalcin Input Parameter: 576211a84d6SLisandro Dalcin + ts - timestepping context 577211a84d6SLisandro Dalcin - use - flag to use adaptivity 578211a84d6SLisandro Dalcin 579211a84d6SLisandro Dalcin Options Database: 580211a84d6SLisandro Dalcin . -ts_bdf_adapt 581211a84d6SLisandro Dalcin 582211a84d6SLisandro Dalcin Level: intermediate 583211a84d6SLisandro Dalcin 584211a84d6SLisandro Dalcin .seealso: TSAdapt 585211a84d6SLisandro Dalcin @*/ 586211a84d6SLisandro Dalcin PetscErrorCode TSBDFUseAdapt(TS ts,PetscBool use) 587211a84d6SLisandro Dalcin { 588211a84d6SLisandro Dalcin PetscErrorCode ierr; 589211a84d6SLisandro Dalcin 590211a84d6SLisandro Dalcin PetscFunctionBegin; 591211a84d6SLisandro Dalcin PetscValidHeaderSpecific(ts,TS_CLASSID,1); 592211a84d6SLisandro Dalcin PetscValidLogicalCollectiveBool(ts,use,2); 593211a84d6SLisandro Dalcin ierr = PetscTryMethod(ts,"TSBDFUseAdapt_C",(TS,PetscBool),(ts,use));CHKERRQ(ierr); 594211a84d6SLisandro Dalcin PetscFunctionReturn(0); 595211a84d6SLisandro Dalcin } 596