1 /* 2 * eimex.c 3 * 4 * Created on: Jun 21, 2012 5 * Written by Hong Zhang (zhang@vt.edu), Virginia Tech 6 * 7 * Contributed by Emil Constantinescu (emconsta@mcs.anl.gov), Argonne National Laboratory 8 */ 9 /* 10 Code for timestepping with Extrapolated IMEX methods 11 12 Notes: 13 The general system is written as 14 15 G(t,X,Xdot) = F(t,X) 16 17 where G represents the stiff part of the physics and F represents the non-stiff part. 18 This method is designed to be linearly implicit on G and can use an approximate and lagged Jacobian. 19 20 Another common form for the system is 21 22 y'=f(x)+g(x) 23 24 The relationship between F,G and f,g is 25 26 G = y'-g(x), F = f(x) 27 28 References 29 E. Constantinescu and A. Sandu, Extrapolated implicit-explicit time stepping, SIAM Journal on Scientific 30 Computing, 31 (2010), pp. 4452–4477. 31 32 */ 33 #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 34 35 static const PetscInt TSEIMEXDefault = 3; 36 37 typedef struct { 38 PetscInt row_ind; /* Return the term T[row_ind][col_ind]*/ 39 PetscInt col_ind; /* Return the term T[row_ind][col_ind]*/ 40 PetscInt nstages; /* Numbers of stages in current scheme*/ 41 PetscInt max_rows; /* Maximum number of rows */ 42 PetscInt *N; /* Harmonic sequence N[max_rows] */ 43 Vec Y; /* States computed during the step, used to complete the step */ 44 Vec Z; /* For shift*(Y-Z) */ 45 Vec *T; /* Working table, size determined by nstages */ 46 Vec YdotRHS; /* f(x) Work vector holding YdotRHS during residual evaluation */ 47 Vec YdotI; /* xdot-g(x) Work vector holding YdotI = G(t,x,xdot) when xdot =0 */ 48 Vec Ydot; /* f(x)+g(x) Work vector*/ 49 Vec VecSolPrev; /* Work vector holding the solution from the previous step (used for interpolation)*/ 50 // PetscScalar *work; /* Scalar work space of length number of stages, used to prepare VecMAXPY() */ 51 PetscReal shift; 52 PetscReal ctime; 53 PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */ 54 PetscBool ord_adapt; /*order adapativity*/ 55 TSStepStatus status; 56 } TS_EIMEX; 57 58 /* This function is pure */ 59 static PetscInt Map(PetscInt i, PetscInt j, PetscInt s) 60 { 61 return ((2*s-j+1)*j/2+i-j); 62 } 63 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "TSEvaluateStep_EIMEX" 67 static PetscErrorCode TSEvaluateStep_EIMEX(TS ts,PetscInt order,Vec X,PetscBool *done) 68 { 69 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 70 const PetscInt ns = ext->nstages; 71 PetscErrorCode ierr; 72 PetscFunctionBegin; 73 ierr = VecCopy(ext->T[Map(ext->row_ind,ext->col_ind,ns)],X);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "TSStage_EIMEX" 80 static PetscErrorCode TSStage_EIMEX(TS ts,PetscInt istage) 81 { 82 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 83 PetscReal h; 84 Vec Y=ext->Y, Z=ext->Z; 85 SNES snes; 86 TSAdapt adapt; 87 PetscInt i,its,lits; 88 PetscBool accept; 89 PetscErrorCode ierr; 90 91 PetscFunctionBegin; 92 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 93 h = ts->time_step/ext->N[istage];/*step size for the istage-th stage*/ 94 ext->shift = 1./h; 95 ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */ 96 ierr = VecCopy(ext->VecSolPrev,Y);CHKERRQ(ierr); /*Take the previous solution as intial step*/ 97 98 for(i=0; i<ext->N[istage]; i++){ 99 ext->ctime = ts->ptime + h*i; 100 ierr = VecCopy(Y,Z);CHKERRQ(ierr);/*Save the solution of the previous substep*/ 101 ierr = SNESSolve(snes,PETSC_NULL,Y);CHKERRQ(ierr); 102 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 103 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 104 ts->snes_its += its; ts->ksp_its += lits; 105 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 106 ierr = TSAdaptCheckStage(adapt,ts,&accept);CHKERRQ(ierr); 107 } 108 109 PetscFunctionReturn(0); 110 } 111 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "TSStep_EIMEX" 115 static PetscErrorCode TSStep_EIMEX(TS ts) 116 { 117 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 118 const PetscInt ns = ext->nstages; 119 Vec *T=ext->T, Y=ext->Y; 120 121 SNES snes; 122 PetscInt i,j; 123 PetscBool accept = PETSC_FALSE; 124 PetscErrorCode ierr; 125 PetscReal alpha,local_error; 126 PetscFunctionBegin; 127 128 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 129 ierr = SNESSetType(snes,"ksponly"); CHKERRQ(ierr); 130 ext->status = TS_STEP_INCOMPLETE; 131 132 ierr = VecCopy(ts->vec_sol,ext->VecSolPrev);CHKERRQ(ierr); 133 134 /* Apply n_j steps of the base method to obtain solutions of T(j,1),1<=j<=s*/ 135 for(j=0; j<ns; j++){ 136 ierr = TSStage_EIMEX(ts,j);CHKERRQ(ierr); 137 ierr = VecCopy(Y,T[j]); CHKERRQ(ierr); 138 } 139 140 for(i=1;i<ns;i++){ 141 for(j=i;j<ns;j++){ 142 alpha = -(PetscReal)ext->N[j]/ext->N[j-i]; 143 ierr = VecAXPBYPCZ(T[Map(j,i,ns)],alpha,1.0,0,T[Map(j,i-1,ns)], 144 T[Map(j-1,i-1,ns)]);/*T[j][i]=alpha*T[j][i-1]+T[j-1][i-1]*/ 145 alpha = 1.0/(1.0 + alpha); 146 ierr = VecScale(T[Map(j,i,ns)],alpha);CHKERRQ(ierr); 147 } 148 } 149 150 ierr = TSEvaluateStep(ts,ns,ts->vec_sol,PETSC_NULL);CHKERRQ(ierr);/*update ts solution */ 151 152 if(ext->ord_adapt && ext->nstages < ext->max_rows){ 153 accept = PETSC_FALSE; 154 while(!accept && ext->nstages < ext->max_rows){ 155 ierr = TSErrorNormWRMS(ts,T[Map(ext->nstages-1,ext->nstages-2,ext->nstages)],&local_error);CHKERRQ(ierr); 156 accept = (local_error < 1.0)? PETSC_TRUE : PETSC_FALSE; 157 158 if(!accept){/* add one more stage*/ 159 ierr = TSStage_EIMEX(ts,ext->nstages);CHKERRQ(ierr); 160 ext->nstages++; ext->row_ind++; ext->col_ind++; 161 /*T table need to be recycled*/ 162 ierr = VecDuplicateVecs(ts->vec_sol,(1+ext->nstages)*ext->nstages/2,&ext->T);CHKERRQ(ierr); 163 for(i=0; i<ext->nstages-1; i++){ 164 for(j=0; j<=i; j++){ 165 ierr = VecCopy(T[Map(i,j,ext->nstages-1)],ext->T[Map(i,j,ext->nstages)]);CHKERRQ(ierr); 166 } 167 } 168 ierr = VecDestroyVecs(ext->nstages*(ext->nstages-1)/2,&T);CHKERRQ(ierr); 169 T = ext->T; /*reset the pointer*/ 170 /*recycling finished, store the new solution*/ 171 ierr = VecCopy(Y,T[ext->nstages-1]); CHKERRQ(ierr); 172 /*extrapolation for the newly added stage*/ 173 for(i=1;i<ext->nstages;i++){ 174 alpha = -(PetscReal)ext->N[ext->nstages-1]/ext->N[ext->nstages-1-i]; 175 ierr = VecAXPBYPCZ(T[Map(ext->nstages-1,i,ext->nstages)],alpha,1.0,0,T[Map(ext->nstages-1,i-1,ext->nstages)],T[Map(ext->nstages-1-1,i-1,ext->nstages)]);/*T[ext->nstages-1][i]=alpha*T[ext->nstages-1][i-1]+T[ext->nstages-1-1][i-1]*/ 176 alpha = 1.0/(1.0 + alpha); 177 ierr = VecScale(T[Map(ext->nstages-1,i,ext->nstages)],alpha);CHKERRQ(ierr); 178 } 179 /*update ts solution */ 180 ierr = TSEvaluateStep(ts,ext->nstages,ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 181 }/*end if !accept*/ 182 }/*end while*/ 183 184 if(ext->nstages == ext->max_rows){ 185 ierr = PetscInfo(ts,"Max number of rows has been used\n");CHKERRQ(ierr); 186 } 187 }/*end if ext->ord_adapt*/ 188 189 ts->ptime += ts->time_step; 190 ts->steps++; 191 ext->status = TS_STEP_COMPLETE; 192 193 if (ext->status != TS_STEP_COMPLETE && !ts->reason) ts->reason = TS_DIVERGED_STEP_REJECTED; 194 PetscFunctionReturn(0); 195 } 196 197 198 /*cubic Hermit spline*/ 199 #undef __FUNCT__ 200 #define __FUNCT__ "TSInterpolate_EIMEX" 201 static PetscErrorCode TSInterpolate_EIMEX(TS ts,PetscReal itime,Vec X) 202 { 203 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 204 PetscReal t,a,b; 205 Vec Y0=ext->VecSolPrev,Y1=ext->Y, 206 Ydot=ext->Ydot,YdotI=ext->YdotI; 207 const PetscReal h = ts->time_step_prev; 208 PetscErrorCode ierr; 209 PetscScalar *x; 210 PetscFunctionBegin; 211 t = (itime -ts->ptime + h)/h; 212 /*YdotI = -f(x)-g(x)*/ 213 214 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); 215 ierr = TSComputeIFunction(ts,ts->ptime-h,Y0,Ydot,YdotI,PETSC_FALSE);CHKERRQ(ierr); 216 217 a = 2.0*pow(t,3) - 3.0*t*t + 1.0; 218 b = -(pow(t,3) - 2.0*t*t + t)*h; 219 ierr = VecAXPBYPCZ(X,a,b,0.0,Y0,YdotI);CHKERRQ(ierr); 220 221 ierr = TSComputeIFunction(ts,ts->ptime,Y1,Ydot,YdotI,PETSC_FALSE);CHKERRQ(ierr); 222 a = -2.0*pow(t,3)+3.0*t*t; 223 b = -(pow(t,3) - t*t)*h; 224 ierr = VecAXPBYPCZ(X,a,b,1.0,Y1,YdotI);CHKERRQ(ierr); 225 226 PetscFunctionReturn(0); 227 } 228 229 230 #undef __FUNCT__ 231 #define __FUNCT__ "TSReset_EIMEX" 232 static PetscErrorCode TSReset_EIMEX(TS ts) 233 { 234 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 235 PetscInt ns; 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 ns = ext->nstages; 240 ierr = VecDestroyVecs((1+ns)*ns/2,&ext->T);CHKERRQ(ierr); 241 ierr = VecDestroy(&ext->Y);CHKERRQ(ierr); 242 ierr = VecDestroy(&ext->Z);CHKERRQ(ierr); 243 ierr = VecDestroy(&ext->YdotRHS);CHKERRQ(ierr); 244 ierr = VecDestroy(&ext->YdotI);CHKERRQ(ierr); 245 ierr = VecDestroy(&ext->Ydot);CHKERRQ(ierr); 246 ierr = VecDestroy(&ext->VecSolPrev);CHKERRQ(ierr); 247 ierr = PetscFree(ext->N);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250 251 #undef __FUNCT__ 252 #define __FUNCT__ "TSDestroy_EIMEX" 253 static PetscErrorCode TSDestroy_EIMEX(TS ts) 254 { 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 ierr = TSReset_EIMEX(ts);CHKERRQ(ierr); 259 ierr = PetscFree(ts->data);CHKERRQ(ierr); 260 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSEIMEXSetMaxRows_C","",PETSC_NULL);CHKERRQ(ierr); 261 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSEIMEXSetRowCol_C","",PETSC_NULL);CHKERRQ(ierr); 262 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSEIMEXSetOrdAdapt_C","",PETSC_NULL);CHKERRQ(ierr); 263 264 PetscFunctionReturn(0); 265 } 266 267 268 #undef __FUNCT__ 269 #define __FUNCT__ "TSEIMEXGetVecs" 270 static PetscErrorCode TSEIMEXGetVecs(TS ts,DM dm,Vec *Z,Vec *Ydot,Vec *YdotI, Vec *YdotRHS) 271 { 272 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 273 PetscErrorCode ierr; 274 275 PetscFunctionBegin; 276 if (Z) { 277 if (dm && dm != ts->dm) { 278 ierr = DMGetNamedGlobalVector(dm,"TSEIMEX_Z",Z);CHKERRQ(ierr); 279 } else *Z = ext->Z; 280 } 281 if (Ydot) { 282 if (dm && dm != ts->dm) { 283 ierr = DMGetNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);CHKERRQ(ierr); 284 } else *Ydot = ext->Ydot; 285 } 286 if (YdotI) { 287 if (dm && dm != ts->dm) { 288 ierr = DMGetNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);CHKERRQ(ierr); 289 } else *YdotI = ext->YdotI; 290 } 291 if (YdotRHS) { 292 if (dm && dm != ts->dm) { 293 ierr = DMGetNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);CHKERRQ(ierr); 294 } else *YdotRHS = ext->YdotRHS; 295 } 296 PetscFunctionReturn(0); 297 } 298 299 300 #undef __FUNCT__ 301 #define __FUNCT__ "TSEIMEXRestoreVecs" 302 static PetscErrorCode TSEIMEXRestoreVecs(TS ts,DM dm,Vec *Z,Vec *Ydot,Vec *YdotI,Vec *YdotRHS) 303 { 304 PetscErrorCode ierr; 305 306 PetscFunctionBegin; 307 if (Z) { 308 if (dm && dm != ts->dm) { 309 ierr = DMRestoreNamedGlobalVector(dm,"TSEIMEX_Z",Z);CHKERRQ(ierr); 310 } 311 } 312 if (Ydot) { 313 if (dm && dm != ts->dm) { 314 ierr = DMRestoreNamedGlobalVector(dm,"TSEIMEX_Ydot",Ydot);CHKERRQ(ierr); 315 } 316 } 317 if (YdotI) { 318 if (dm && dm != ts->dm) { 319 ierr = DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotI",YdotI);CHKERRQ(ierr); 320 } 321 } 322 if (YdotRHS) { 323 if (dm && dm != ts->dm) { 324 ierr = DMRestoreNamedGlobalVector(dm,"TSEIMEX_YdotRHS",YdotRHS);CHKERRQ(ierr); 325 } 326 } 327 PetscFunctionReturn(0); 328 } 329 330 331 /* 332 This defines the nonlinear equation that is to be solved with SNES 333 Fn[t0+Theta*dt, U, (U-U0)*shift] = 0 334 In the case of Backward Euler, Fn = (U-U0)/h-g(t1,U)) 335 Since FormIFunction calculates G = ydot - g(t,y), ydot will be set to (U-U0)/h 336 */ 337 #undef __FUNCT__ 338 #define __FUNCT__ "SNESTSFormFunction_EIMEX" 339 static PetscErrorCode SNESTSFormFunction_EIMEX(SNES snes,Vec X,Vec G,TS ts) 340 { 341 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 342 PetscErrorCode ierr; 343 Vec Ydot,Z; 344 DM dm,dmsave; 345 346 PetscFunctionBegin; 347 ierr = VecZeroEntries(G);CHKERRQ(ierr); 348 349 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 350 ierr = TSEIMEXGetVecs(ts,dm,&Z,&Ydot,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 351 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); 352 dmsave = ts->dm; 353 ts->dm = dm; 354 ierr = TSComputeIFunction(ts,ext->ctime,X,Ydot,G,PETSC_FALSE);CHKERRQ(ierr); 355 /* PETSC_FALSE indicates non-imex, adding explicit RHS to the implicit I function. */ 356 ierr = VecCopy(G,Ydot);CHKERRQ(ierr); 357 ts->dm = dmsave; 358 ierr = TSEIMEXRestoreVecs(ts,dm,&Z,&Ydot,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 359 360 PetscFunctionReturn(0); 361 } 362 363 /* 364 This defined the Jacobian matrix for SNES. Jn = (I/h-g'(t,y)) 365 */ 366 #undef __FUNCT__ 367 #define __FUNCT__ "SNESTSFormJacobian_EIMEX" 368 static PetscErrorCode SNESTSFormJacobian_EIMEX(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 369 { 370 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 371 Vec Ydot; 372 PetscErrorCode ierr; 373 DM dm,dmsave; 374 PetscFunctionBegin; 375 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 376 ierr = TSEIMEXGetVecs(ts,dm,PETSC_NULL,&Ydot,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 377 /* ierr = VecZeroEntries(Ydot);CHKERRQ(ierr);*/ 378 /* ext->Ydot have already been computed in SNESTSFormFunction_EIMEX (SNES guarantees this) */ 379 dmsave = ts->dm; 380 ts->dm = dm; 381 ierr = TSComputeIJacobian(ts,ts->ptime,X,Ydot,ext->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 382 ts->dm = dmsave; 383 ierr = TSEIMEXRestoreVecs(ts,dm,PETSC_NULL,&Ydot,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 384 PetscFunctionReturn(0); 385 } 386 387 #undef __FUNCT__ 388 #define __FUNCT__ "DMCoarsenHook_TSEIMEX" 389 static PetscErrorCode DMCoarsenHook_TSEIMEX(DM fine,DM coarse,void *ctx) 390 { 391 392 PetscFunctionBegin; 393 PetscFunctionReturn(0); 394 } 395 396 #undef __FUNCT__ 397 #define __FUNCT__ "DMRestrictHook_TSEIMEX" 398 static PetscErrorCode DMRestrictHook_TSEIMEX(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse,void *ctx) 399 { 400 TS ts = (TS)ctx; 401 PetscErrorCode ierr; 402 Vec Z,Z_c; 403 404 PetscFunctionBegin; 405 ierr = TSEIMEXGetVecs(ts,fine,&Z,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 406 ierr = TSEIMEXGetVecs(ts,coarse,&Z_c,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 407 ierr = MatRestrict(restrct,Z,Z_c);CHKERRQ(ierr); 408 ierr = VecPointwiseMult(Z_c,rscale,Z_c);CHKERRQ(ierr); 409 ierr = TSEIMEXRestoreVecs(ts,fine,&Z,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 410 ierr = TSEIMEXRestoreVecs(ts,coarse,&Z_c,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 411 PetscFunctionReturn(0); 412 } 413 414 415 #undef __FUNCT__ 416 #define __FUNCT__ "TSSetUp_EIMEX" 417 static PetscErrorCode TSSetUp_EIMEX(TS ts) 418 { 419 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 420 PetscErrorCode ierr; 421 DM dm; 422 PetscFunctionBegin; 423 424 if (!ext->N){ /*ext->max_rows not set*/ 425 ierr = TSEIMEXSetMaxRows(ts,TSEIMEXDefault);CHKERRQ(ierr); 426 } 427 if(-1 == ext->row_ind && -1 == ext->col_ind){ 428 ierr = TSEIMEXSetRowCol(ts,ext->max_rows,ext->max_rows);CHKERRQ(ierr); 429 }else{/*ext->row_ind and col_ind already set*/ 430 if(ext->ord_adapt){ 431 ierr = PetscInfo(ts,"Order adaptivity is enabled and TSEIMEXSetRowCol or -ts_eimex_row_col option will take no effect\n");CHKERRQ(ierr); 432 } 433 } 434 435 if(ext->ord_adapt){ 436 ext->nstages = 2; /*Start with the 2-stage scheme*/ 437 ierr = TSEIMEXSetRowCol(ts,ext->nstages,ext->nstages);CHKERRQ(ierr); 438 }else{ 439 ext->nstages = ext->max_rows; /*by default nstages is the same as max_rows, this can be changed by setting order adaptivity */ 440 } 441 442 ierr = VecDuplicateVecs(ts->vec_sol,(1+ext->nstages)*ext->nstages/2,&ext->T);CHKERRQ(ierr);/*full T table*/ 443 ierr = VecDuplicate(ts->vec_sol,&ext->YdotI);CHKERRQ(ierr); 444 ierr = VecDuplicate(ts->vec_sol,&ext->YdotRHS);CHKERRQ(ierr); 445 ierr = VecDuplicate(ts->vec_sol,&ext->Ydot);CHKERRQ(ierr); 446 ierr = VecDuplicate(ts->vec_sol,&ext->VecSolPrev);CHKERRQ(ierr); 447 ierr = VecDuplicate(ts->vec_sol,&ext->Y);CHKERRQ(ierr); 448 ierr = VecDuplicate(ts->vec_sol,&ext->Z);CHKERRQ(ierr); 449 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 450 if (dm) { 451 ierr = DMCoarsenHookAdd(dm,DMCoarsenHook_TSEIMEX,DMRestrictHook_TSEIMEX,ts);CHKERRQ(ierr); 452 } 453 PetscFunctionReturn(0); 454 } 455 456 #undef __FUNCT__ 457 #define __FUNCT__ "TSSetFromOptions_EIMEX" 458 static PetscErrorCode TSSetFromOptions_EIMEX(TS ts) 459 { 460 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 461 PetscErrorCode ierr; 462 PetscInt tindex[2]={TSEIMEXDefault,TSEIMEXDefault}; 463 PetscInt np = 2, nrows=TSEIMEXDefault; 464 PetscFunctionBegin; 465 ierr = PetscOptionsHead("EIMEX ODE solver options");CHKERRQ(ierr); 466 { 467 PetscBool flg; 468 flg = PETSC_FALSE; 469 ierr = PetscOptionsInt("-ts_eimex_max_rows","Define the maximum number of rows used","TSEIMEXSetMaxRows",nrows,&nrows,&flg);CHKERRQ(ierr); /*default value 3*/ 470 if(flg){ 471 ierr = TSEIMEXSetMaxRows(ts,nrows);CHKERRQ(ierr); 472 } 473 ierr = PetscOptionsIntArray("-ts_eimex_row_col","Return the specific term in the T table","TSEIMEXSetRowCol",tindex,&np,&flg);CHKERRQ(ierr); 474 if(flg){ 475 ierr = TSEIMEXSetRowCol(ts,tindex[0],tindex[1]);CHKERRQ(ierr); 476 } 477 ierr = PetscOptionsBool("-ts_eimex_order_adapt","Solve the problem with adaptive order","TSEIMEXSetOrdAdapt",ext->ord_adapt,&ext->ord_adapt,PETSC_NULL);CHKERRQ(ierr); 478 ierr = SNESSetFromOptions(ts->snes);CHKERRQ(ierr); 479 480 /*prevent users from changing SNESType*/ 481 } 482 ierr = PetscOptionsTail();CHKERRQ(ierr); 483 PetscFunctionReturn(0); 484 } 485 486 #undef __FUNCT__ 487 #define __FUNCT__ "TSView_EIMEX" 488 static PetscErrorCode TSView_EIMEX(TS ts,PetscViewer viewer) 489 { 490 /* TS_EIMEX *ext = (TS_EIMEX*)ts->data;*/ 491 PetscBool iascii; 492 PetscErrorCode ierr; 493 494 PetscFunctionBegin; 495 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 496 if (iascii) { 497 ierr = PetscViewerASCIIPrintf(viewer," EIMEX\n");CHKERRQ(ierr); 498 } 499 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 500 PetscFunctionReturn(0); 501 } 502 503 504 #undef __FUNCT__ 505 #define __FUNCT__ "TSEIMEXSetMaxRows" 506 /*@C 507 TSEIMEXSetMaxRows - Set the maximum number of rows for EIMEX schemes 508 509 Logically collective 510 511 Input Parameter: 512 + ts - timestepping context 513 - nrows - maximum number of rows 514 515 Level: intermediate 516 517 .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX 518 @*/ 519 PetscErrorCode TSEIMEXSetMaxRows(TS ts, PetscInt nrows) 520 { 521 PetscErrorCode ierr; 522 PetscFunctionBegin; 523 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 524 ierr = PetscTryMethod(ts,"TSEIMEXSetMaxRows_C",(TS,PetscInt),(ts,nrows));CHKERRQ(ierr); 525 PetscFunctionReturn(0); 526 } 527 528 529 #undef __FUNCT__ 530 #define __FUNCT__ "TSEIMEXSetRowCol" 531 /*@C 532 TSEIMEXSetRowCol - Set the type index in the T table for the return value 533 534 Logically collective 535 536 Input Parameter: 537 + ts - timestepping context 538 - tindex - index in the T table 539 540 Level: intermediate 541 542 .seealso: TSEIMEXSetMaxRows(), TSEIMEXSetOrdAdapt(), TSEIMEX 543 @*/ 544 PetscErrorCode TSEIMEXSetRowCol(TS ts, PetscInt row, PetscInt col) 545 { 546 PetscErrorCode ierr; 547 PetscFunctionBegin; 548 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 549 ierr = PetscTryMethod(ts,"TSEIMEXSetRowCol_C",(TS,PetscInt, PetscInt),(ts,row,col));CHKERRQ(ierr); 550 PetscFunctionReturn(0); 551 } 552 553 554 #undef __FUNCT__ 555 #define __FUNCT__ "TSEIMEXSetOrdAdapt" 556 /*@C 557 TSEIMEXSetOrdAdapt - Set the order adaptativity 558 559 Logically collective 560 561 Input Parameter: 562 + ts - timestepping context 563 - tindex - index in the T table 564 565 Level: intermediate 566 567 .seealso: TSEIMEXSetRowCol(), TSEIMEXSetOrdAdapt(), TSEIMEX 568 @*/ 569 PetscErrorCode TSEIMEXSetOrdAdapt(TS ts, PetscBool flg) 570 { 571 PetscErrorCode ierr; 572 PetscFunctionBegin; 573 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 574 ierr = PetscTryMethod(ts,"TSEIMEXSetOrdAdapt_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 575 PetscFunctionReturn(0); 576 } 577 578 579 EXTERN_C_BEGIN 580 #undef __FUNCT__ 581 #define __FUNCT__ "TSEIMEXSetMaxRows_EIMEX" 582 PetscErrorCode TSEIMEXSetMaxRows_EIMEX(TS ts,PetscInt nrows) 583 { 584 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 585 PetscErrorCode ierr; 586 PetscInt i; 587 588 PetscFunctionBegin; 589 if (nrows < 0 || nrows > 100) SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Max number of rows (current value %D) should be an integer number between 1 and 100\n",nrows); 590 if(ext->N){ 591 ierr = PetscFree(ext->N);CHKERRQ(ierr); 592 } 593 ext->max_rows = nrows; 594 ierr = PetscMalloc(nrows*sizeof(PetscInt),&ext->N);CHKERRQ(ierr); 595 for(i=0;i<nrows;i++) ext->N[i]=i+1; 596 PetscFunctionReturn(0); 597 } 598 599 #undef __FUNCT__ 600 #define __FUNCT__ "TSEIMEXSetRowCol_EIMEX" 601 PetscErrorCode TSEIMEXSetRowCol_EIMEX(TS ts,PetscInt row,PetscInt col) 602 { 603 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 604 605 PetscFunctionBegin; 606 if (row < 1 || col < 1) SETERRQ2(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"The row or column index (current value %d,%d) should not be less than 1 \n",row,col); 607 if (row > ext->max_rows || col > ext->max_rows) SETERRQ3(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"The row or column index (current value %d,%d) exceeds the maximum number of rows %d\n",row,col,ext->max_rows); 608 if (col > row) SETERRQ2(((PetscObject)ts)->comm,PETSC_ERR_ARG_OUTOFRANGE,"The column index (%d) exceeds the row index (%d)\n",col,row); 609 610 ext->row_ind = row - 1; 611 ext->col_ind = col - 1; /*Array index in C starts from 0*/ 612 PetscFunctionReturn(0); 613 } 614 615 #undef __FUNCT__ 616 #define __FUNCT__ "TSEIMEXSetOrdAdapt_EIMEX" 617 PetscErrorCode TSEIMEXSetOrdAdapt_EIMEX(TS ts,PetscBool flg) 618 { 619 TS_EIMEX *ext = (TS_EIMEX*)ts->data; 620 PetscFunctionBegin; 621 ext->ord_adapt = flg; 622 PetscFunctionReturn(0); 623 } 624 EXTERN_C_END 625 626 /* ------------------------------------------------------------ */ 627 /*MC 628 TSEIMEX - ODE solver using extrapolated IMEX schemes 629 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 630 631 Notes: 632 The default is a 3-stage scheme, it can be changed with TSEIMEXSetMaxRows() or -ts_eimex_max_rows 633 634 This method currently only works with ODE, for which the stiff part G(t,X,Xdot) has the form Xdot + Ghat(t,X). 635 636 Level: beginner 637 638 .seealso: TSCreate(), TS 639 M*/ 640 EXTERN_C_BEGIN 641 #undef __FUNCT__ 642 #define __FUNCT__ "TSCreate_EIMEX" 643 PetscErrorCode TSCreate_EIMEX(TS ts) 644 { 645 TS_EIMEX *ext; 646 PetscErrorCode ierr; 647 648 PetscFunctionBegin; 649 650 ts->ops->reset = TSReset_EIMEX; 651 ts->ops->destroy = TSDestroy_EIMEX; 652 ts->ops->view = TSView_EIMEX; 653 ts->ops->setup = TSSetUp_EIMEX; 654 ts->ops->step = TSStep_EIMEX; 655 ts->ops->interpolate = TSInterpolate_EIMEX; 656 ts->ops->evaluatestep = TSEvaluateStep_EIMEX; 657 ts->ops->setfromoptions = TSSetFromOptions_EIMEX; 658 ts->ops->snesfunction = SNESTSFormFunction_EIMEX; 659 ts->ops->snesjacobian = SNESTSFormJacobian_EIMEX; 660 661 ierr = PetscNewLog(ts,TS_EIMEX,&ext);CHKERRQ(ierr); 662 ts->data = (void*)ext; 663 664 ext->ord_adapt = PETSC_FALSE; /*By default, no order adapativity*/ 665 ext->row_ind = -1; 666 ext->col_ind = -1; 667 ext->max_rows = TSEIMEXDefault; 668 ext->nstages = TSEIMEXDefault; 669 670 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSEIMEXSetMaxRows_C","TSEIMEXSetMaxRows_EIMEX",TSEIMEXSetMaxRows_EIMEX);CHKERRQ(ierr); 671 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSEIMEXSetRowCol_C","TSEIMEXSetRowCol_EIMEX",TSEIMEXSetRowCol_EIMEX);CHKERRQ(ierr); 672 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSEIMEXSetOrdAdapt_C","TSEIMEXSetOrdAdapt_EIMEX",TSEIMEXSetOrdAdapt_EIMEX);CHKERRQ(ierr); 673 PetscFunctionReturn(0); 674 } 675 EXTERN_C_END 676