1 /* 2 Code for timestepping with additive Runge-Kutta IMEX method 3 4 Notes: 5 The general system is written as 6 7 F(t,X,Xdot) = G(t,X) 8 9 where F represents the stiff part of the physics and G represents the non-stiff part. 10 11 */ 12 #include <private/tsimpl.h> /*I "petscts.h" I*/ 13 14 static const TSARKIMEXType TSARKIMEXDefault = TSARKIMEX2D; 15 static PetscBool TSARKIMEXRegisterAllCalled; 16 static PetscBool TSARKIMEXPackageInitialized; 17 18 typedef struct _ARKTableau *ARKTableau; 19 struct _ARKTableau { 20 char *name; 21 PetscInt order; 22 PetscInt s; 23 PetscReal *At,*bt,*ct; 24 PetscReal *A,*b,*c; /* Non-stiff tableau */ 25 }; 26 typedef struct _ARKTableauLink *ARKTableauLink; 27 struct _ARKTableauLink { 28 struct _ARKTableau tab; 29 ARKTableauLink next; 30 }; 31 static ARKTableauLink ARKTableauList; 32 33 typedef struct { 34 ARKTableau tableau; 35 Vec *Y; /* States computed during the step */ 36 Vec *YdotI; /* Time derivatives for the stiff part */ 37 Vec *YdotRHS; /* Function evaluations for the non-stiff part */ 38 Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 39 Vec Work; /* Generic work vector */ 40 Vec Z; /* Ydot = shift(Y-Z) */ 41 PetscScalar *work; /* Scalar work */ 42 PetscReal shift; 43 PetscReal stage_time; 44 } TS_ARKIMEX; 45 46 #undef __FUNCT__ 47 #define __FUNCT__ "TSARKIMEXRegisterAll" 48 /*@C 49 TSARKIMEXRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSARKIMEX 50 51 Not Collective 52 53 Level: advanced 54 55 .keywords: TS, TSARKIMEX, register, all 56 57 .seealso: TSARKIMEXRegisterDestroy() 58 @*/ 59 PetscErrorCode TSARKIMEXRegisterAll(void) 60 { 61 PetscErrorCode ierr; 62 63 PetscFunctionBegin; 64 if (TSARKIMEXRegisterAllCalled) PetscFunctionReturn(0); 65 TSARKIMEXRegisterAllCalled = PETSC_TRUE; 66 67 { 68 const PetscReal 69 A[3][3] = {{0,0,0}, 70 {0.41421356237309504880,0,0}, 71 {0.75,0.25,0}}, 72 At[3][3] = {{0,0,0}, 73 {0.12132034355964257320,0.29289321881345247560,0}, 74 {0.20710678118654752440,0.50000000000000000000,0.29289321881345247560}}; 75 ierr = TSARKIMEXRegister(TSARKIMEX2D,2,3,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 76 } 77 { 78 const PetscReal s2 = sqrt(2), 79 A[3][3] = {{0,0,0}, 80 {2-s2,0,0}, 81 {(3-2*s2)/6,(3+2*s2)/6,0}}, 82 At[3][3] = {{0,0,0}, 83 {1-1/s2,1-1/s2,0}, 84 {1/(2*s2),1/(2*s2),1-1/s2}}; 85 ierr = TSARKIMEXRegister(TSARKIMEX2E,2,3,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 86 } 87 { 88 const PetscReal 89 A[4][4] = {{0,0,0,0}, 90 {1767732205903./2027836641118,0,0,0}, 91 {5535828885825./10492691773637,788022342437./10882634858940,0,0}, 92 {6485989280629./16251701735622,-4246266847089./9704473918619,10755448449292./10357097424841,0}}, 93 At[4][4] = {{0,0,0,0}, 94 {1767732205903./4055673282236,1767732205903./4055673282236,0,0}, 95 {2746238789719./10658868560708,-640167445237./6845629431997,1767732205903./4055673282236,0}, 96 {1471266399579./7840856788654,-4482444167858./7529755066697,11266239266428./11593286722821,1767732205903./4055673282236}}; 97 ierr = TSARKIMEXRegister(TSARKIMEX3,3,4,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 98 } 99 { 100 const PetscReal 101 A[6][6] = {{0,0,0,0,0,0}, 102 {1./2,0,0,0,0,0}, 103 {13861./62500,6889./62500,0,0,0,0}, 104 {-116923316275./2393684061468,-2731218467317./15368042101831,9408046702089./11113171139209,0,0,0}, 105 {-451086348788./2902428689909,-2682348792572./7519795681897,12662868775082./11960479115383,3355817975965./11060851509271,0,0}, 106 {647845179188./3216320057751,73281519250./8382639484533,552539513391./3454668386233,3354512671639./8306763924573,4040./17871,0}}, 107 At[6][6] = {{0,0,0,0,0,0}, 108 {1./4,1./4,0,0,0,0}, 109 {8611./62500,-1743./31250,1./4,0,0,0}, 110 {5012029./34652500,-654441./2922500,174375./388108,1./4,0,0}, 111 {15267082809./155376265600,-71443401./120774400,730878875./902184768,2285395./8070912,1./4,0}, 112 {82889./524892,0,15625./83664,69875./102672,-2260./8211,1./4}}; 113 ierr = TSARKIMEXRegister(TSARKIMEX4,4,6,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 114 } 115 { 116 const PetscReal 117 A[8][8] = {{0,0,0,0,0,0,0,0}, 118 {41./100,0,0,0,0,0,0,0}, 119 {367902744464./2072280473677,677623207551./8224143866563,0,0,0,0,0,0}, 120 {1268023523408./10340822734521,0,1029933939417./13636558850479,0,0,0,0,0}, 121 {14463281900351./6315353703477,0,66114435211212./5879490589093,-54053170152839./4284798021562,0,0,0,0}, 122 {14090043504691./34967701212078,0,15191511035443./11219624916014,-18461159152457./12425892160975,-281667163811./9011619295870,0,0,0}, 123 {19230459214898./13134317526959,0,21275331358303./2942455364971,-38145345988419./4862620318723,-1./8,-1./8,0,0}, 124 {-19977161125411./11928030595625,0,-40795976796054./6384907823539,177454434618887./12078138498510,782672205425./8267701900261,-69563011059811./9646580694205,7356628210526./4942186776405,0}}, 125 At[8][8] = {{0,0,0,0,0,0,0,0}, 126 {41./200,41./200,0,0,0,0,0,0}, 127 {41./400,-567603406766./11931857230679,41./200,0,0,0,0,0}, 128 {683785636431./9252920307686,0,-110385047103./1367015193373,41./200,0,0,0,0}, 129 {3016520224154./10081342136671,0,30586259806659./12414158314087,-22760509404356./11113319521817,41./200,0,0,0}, 130 {218866479029./1489978393911,0,638256894668./5436446318841,-1179710474555./5321154724896,-60928119172./8023461067671,41./200,0,0}, 131 {1020004230633./5715676835656,0,25762820946817./25263940353407,-2161375909145./9755907335909,-211217309593./5846859502534,-4269925059573./7827059040749,41./200,0}, 132 {-872700587467./9133579230613,0,0,22348218063261./9555858737531,-1143369518992./8141816002931,-39379526789629./19018526304540,32727382324388./42900044865799,41./200}}; 133 ierr = TSARKIMEXRegister(TSARKIMEX5,5,8,&At[0][0],PETSC_NULL,PETSC_NULL,&A[0][0],PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 134 } 135 136 PetscFunctionReturn(0); 137 } 138 139 #undef __FUNCT__ 140 #define __FUNCT__ "TSARKIMEXRegisterDestroy" 141 /*@C 142 TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by TSARKIMEXRegister(). 143 144 Not Collective 145 146 Level: advanced 147 148 .keywords: TSARKIMEX, register, destroy 149 .seealso: TSARKIMEXRegister(), TSARKIMEXRegisterAll(), TSARKIMEXRegisterDynamic() 150 @*/ 151 PetscErrorCode TSARKIMEXRegisterDestroy(void) 152 { 153 PetscErrorCode ierr; 154 ARKTableauLink link; 155 156 PetscFunctionBegin; 157 while ((link = ARKTableauList)) { 158 ARKTableau t = &link->tab; 159 ARKTableauList = link->next; 160 ierr = PetscFree6(t->At,t->bt,t->ct,t->A,t->b,t->c);CHKERRQ(ierr); 161 ierr = PetscFree(t->name);CHKERRQ(ierr); 162 ierr = PetscFree(link);CHKERRQ(ierr); 163 } 164 TSARKIMEXRegisterAllCalled = PETSC_FALSE; 165 PetscFunctionReturn(0); 166 } 167 168 #undef __FUNCT__ 169 #define __FUNCT__ "TSARKIMEXInitializePackage" 170 /*@C 171 TSARKIMEXInitializePackage - This function initializes everything in the TSARKIMEX package. It is called 172 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_ARKIMEX() 173 when using static libraries. 174 175 Input Parameter: 176 path - The dynamic library path, or PETSC_NULL 177 178 Level: developer 179 180 .keywords: TS, TSARKIMEX, initialize, package 181 .seealso: PetscInitialize() 182 @*/ 183 PetscErrorCode TSARKIMEXInitializePackage(const char path[]) 184 { 185 PetscErrorCode ierr; 186 187 PetscFunctionBegin; 188 if (TSARKIMEXPackageInitialized) PetscFunctionReturn(0); 189 TSARKIMEXPackageInitialized = PETSC_TRUE; 190 ierr = TSARKIMEXRegisterAll();CHKERRQ(ierr); 191 ierr = PetscRegisterFinalize(TSARKIMEXFinalizePackage);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 195 #undef __FUNCT__ 196 #define __FUNCT__ "TSARKIMEXFinalizePackage" 197 /*@C 198 TSARKIMEXFinalizePackage - This function destroys everything in the TSARKIMEX package. It is 199 called from PetscFinalize(). 200 201 Level: developer 202 203 .keywords: Petsc, destroy, package 204 .seealso: PetscFinalize() 205 @*/ 206 PetscErrorCode TSARKIMEXFinalizePackage(void) 207 { 208 PetscErrorCode ierr; 209 210 PetscFunctionBegin; 211 TSARKIMEXPackageInitialized = PETSC_FALSE; 212 ierr = TSARKIMEXRegisterDestroy();CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 #undef __FUNCT__ 217 #define __FUNCT__ "TSARKIMEXRegister" 218 PetscErrorCode TSARKIMEXRegister(const TSARKIMEXType name,PetscInt order,PetscInt s, 219 const PetscReal At[],const PetscReal bt[],const PetscReal ct[], 220 const PetscReal A[],const PetscReal b[],const PetscReal c[]) 221 { 222 PetscErrorCode ierr; 223 ARKTableauLink link; 224 ARKTableau t; 225 PetscInt i,j; 226 227 PetscFunctionBegin; 228 ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 229 t = &link->tab; 230 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 231 t->order = order; 232 t->s = s; 233 ierr = PetscMalloc6(s*s,PetscReal,&t->At,s,PetscReal,&t->bt,s,PetscReal,&t->ct,s*s,PetscReal,&t->A,s,PetscReal,&t->b,s,PetscReal,&t->c);CHKERRQ(ierr); 234 ierr = PetscMemcpy(t->At,At,s*s*sizeof(At[0]));CHKERRQ(ierr); 235 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 236 if (bt) {ierr = PetscMemcpy(t->bt,bt,s*sizeof(bt[0]));CHKERRQ(ierr);} 237 else for (i=0; i<s; i++) t->bt[i] = At[(s-1)*s+i]; 238 if (b) {ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr);} 239 else for (i=0; i<s; i++) t->b[i] = At[(s-1)*s+i]; 240 if (ct) {ierr = PetscMemcpy(t->ct,ct,s*sizeof(ct[0]));CHKERRQ(ierr);} 241 else for (i=0; i<s; i++) for (j=0,t->ct[i]=0; j<s; j++) t->ct[i] += At[i*s+j]; 242 if (c) {ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr);} 243 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 244 link->next = ARKTableauList; 245 ARKTableauList = link; 246 PetscFunctionReturn(0); 247 } 248 249 #undef __FUNCT__ 250 #define __FUNCT__ "TSStep_ARKIMEX" 251 static PetscErrorCode TSStep_ARKIMEX(TS ts) 252 { 253 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 254 ARKTableau tab = ark->tableau; 255 const PetscInt s = tab->s; 256 const PetscReal *At = tab->At,*A = tab->A,*bt = tab->bt,*b = tab->b,*ct = tab->ct,*c = tab->c; 257 PetscReal *w = ark->work; 258 Vec *Y = ark->Y,*YdotI = ark->YdotI,*YdotRHS = ark->YdotRHS,Ydot = ark->Ydot,W = ark->Work,Z = ark->Z; 259 SNES snes; 260 PetscInt i,j,its,lits; 261 PetscReal h,t; 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 266 h = ts->time_step = ts->next_time_step; 267 t = ts->ptime; 268 269 for (i=0; i<s; i++) { 270 if (At[i*s+i] == 0) { /* This stage is explicit */ 271 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 272 for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 273 ierr = VecMAXPY(Y[i],i,w,YdotI);CHKERRQ(ierr); 274 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 275 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 276 } else { 277 ark->stage_time = t + h*ct[i]; 278 ark->shift = 1./(h*At[i*s+i]); 279 /* Affine part */ 280 ierr = VecZeroEntries(W);CHKERRQ(ierr); 281 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 282 ierr = VecMAXPY(W,i,w,YdotRHS);CHKERRQ(ierr); 283 /* Ydot = shift*(Y-Z) */ 284 ierr = VecCopy(ts->vec_sol,Z);CHKERRQ(ierr); 285 for (j=0; j<i; j++) w[j] = h*At[i*s+j]; 286 ierr = VecMAXPY(Z,i,w,YdotRHS);CHKERRQ(ierr); 287 /* Initial guess taken from last stage */ 288 ierr = VecCopy(i>0?Y[i-1]:ts->vec_sol,Y[i]);CHKERRQ(ierr); 289 ierr = SNESSolve(snes,W,Y[i]);CHKERRQ(ierr); 290 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 291 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 292 ts->nonlinear_its += its; ts->linear_its += lits; 293 } 294 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); 295 ierr = TSComputeIFunction(ts,t+h*ct[i],Y[i],Ydot,YdotI[i],PETSC_TRUE);CHKERRQ(ierr); 296 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 297 } 298 for (j=0; j<s; j++) w[j] = -h*bt[j]; 299 ierr = VecMAXPY(ts->vec_sol,s,w,YdotI);CHKERRQ(ierr); 300 for (j=0; j<s; j++) w[j] = h*b[j]; 301 ierr = VecMAXPY(ts->vec_sol,s,w,YdotRHS);CHKERRQ(ierr); 302 303 ts->ptime += ts->time_step; 304 ts->next_time_step = ts->time_step; 305 ts->steps++; 306 PetscFunctionReturn(0); 307 } 308 309 /*------------------------------------------------------------*/ 310 #undef __FUNCT__ 311 #define __FUNCT__ "TSReset_ARKIMEX" 312 static PetscErrorCode TSReset_ARKIMEX(TS ts) 313 { 314 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 315 PetscInt s; 316 PetscErrorCode ierr; 317 318 PetscFunctionBegin; 319 if (!ark->tableau) PetscFunctionReturn(0); 320 s = ark->tableau->s; 321 ierr = VecDestroyVecs(s,&ark->Y);CHKERRQ(ierr); 322 ierr = VecDestroyVecs(s,&ark->YdotI);CHKERRQ(ierr); 323 ierr = VecDestroyVecs(s,&ark->YdotRHS);CHKERRQ(ierr); 324 ierr = VecDestroy(&ark->Ydot);CHKERRQ(ierr); 325 ierr = VecDestroy(&ark->Work);CHKERRQ(ierr); 326 ierr = VecDestroy(&ark->Z);CHKERRQ(ierr); 327 ierr = PetscFree(ark->work);CHKERRQ(ierr); 328 PetscFunctionReturn(0); 329 } 330 331 #undef __FUNCT__ 332 #define __FUNCT__ "TSDestroy_ARKIMEX" 333 static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 334 { 335 PetscErrorCode ierr; 336 337 PetscFunctionBegin; 338 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 339 ierr = PetscFree(ts->data);CHKERRQ(ierr); 340 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","",PETSC_NULL);CHKERRQ(ierr); 341 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","",PETSC_NULL);CHKERRQ(ierr); 342 PetscFunctionReturn(0); 343 } 344 345 /* 346 This defines the nonlinear equation that is to be solved with SNES 347 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 348 */ 349 #undef __FUNCT__ 350 #define __FUNCT__ "SNESTSFormFunction_ARKIMEX" 351 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts) 352 { 353 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 354 PetscErrorCode ierr; 355 356 PetscFunctionBegin; 357 ierr = VecAXPBYPCZ(ark->Ydot,-ark->shift,ark->shift,0,ark->Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */ 358 ierr = TSComputeIFunction(ts,ark->stage_time,X,ark->Ydot,F,PETSC_TRUE);CHKERRQ(ierr); 359 PetscFunctionReturn(0); 360 } 361 362 #undef __FUNCT__ 363 #define __FUNCT__ "SNESTSFormJacobian_ARKIMEX" 364 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 365 { 366 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 367 PetscErrorCode ierr; 368 369 PetscFunctionBegin; 370 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 371 ierr = TSComputeIJacobian(ts,ark->stage_time,X,ark->Ydot,ark->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 372 PetscFunctionReturn(0); 373 } 374 375 #undef __FUNCT__ 376 #define __FUNCT__ "TSSetUp_ARKIMEX" 377 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 378 { 379 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 380 ARKTableau tab = ark->tableau; 381 PetscInt s = tab->s; 382 PetscErrorCode ierr; 383 384 PetscFunctionBegin; 385 if (!ark->tableau) { 386 ierr = TSSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 387 } 388 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr); 389 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr); 390 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr); 391 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr); 392 ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr); 393 ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr); 394 ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr); 395 PetscFunctionReturn(0); 396 } 397 /*------------------------------------------------------------*/ 398 399 #undef __FUNCT__ 400 #define __FUNCT__ "TSSetFromOptions_ARKIMEX" 401 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts) 402 { 403 PetscErrorCode ierr; 404 char arktype[256]; 405 406 PetscFunctionBegin; 407 ierr = PetscOptionsHead("ARKIMEX ODE solver options");CHKERRQ(ierr); 408 { 409 ARKTableauLink link; 410 PetscInt count,choice; 411 PetscBool flg; 412 const char **namelist; 413 ierr = PetscStrncpy(arktype,TSARKIMEXDefault,sizeof arktype);CHKERRQ(ierr); 414 for (link=ARKTableauList,count=0; link; link=link->next,count++) ; 415 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 416 for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 417 ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 418 ierr = TSARKIMEXSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 419 ierr = PetscFree(namelist);CHKERRQ(ierr); 420 } 421 ierr = PetscOptionsTail();CHKERRQ(ierr); 422 PetscFunctionReturn(0); 423 } 424 425 #undef __FUNCT__ 426 #define __FUNCT__ "PetscFormatRealArray" 427 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 428 { 429 int i,left,count; 430 char *p; 431 432 PetscFunctionBegin; 433 for (i=0,p=buf,left=(int)len; i<n; i++) { 434 count = snprintf(p,left,fmt,x[i]); 435 if (count < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"snprintf()"); 436 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 437 left -= count; 438 p += count; 439 *p++ = ' '; 440 } 441 p[i ? 0 : -1] = 0; 442 PetscFunctionReturn(0); 443 } 444 445 #undef __FUNCT__ 446 #define __FUNCT__ "TSView_ARKIMEX" 447 static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer) 448 { 449 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 450 ARKTableau tab = ark->tableau; 451 PetscBool iascii; 452 PetscErrorCode ierr; 453 454 PetscFunctionBegin; 455 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 456 if (iascii) { 457 const TSARKIMEXType arktype; 458 char buf[512]; 459 ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr); 460 ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 461 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ct);CHKERRQ(ierr); 462 ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa ct = %s\n",buf);CHKERRQ(ierr); 463 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->c);CHKERRQ(ierr); 464 ierr = PetscViewerASCIIPrintf(viewer," Nonstiff abscissa c = %s\n",buf);CHKERRQ(ierr); 465 } else { 466 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for TS_ARKIMEX",((PetscObject)viewer)->type_name); 467 } 468 PetscFunctionReturn(0); 469 } 470 471 #undef __FUNCT__ 472 #define __FUNCT__ "TSARKIMEXSetType" 473 /*@C 474 TSARKIMEXSetType - Set the type of ARK IMEX scheme 475 476 Logically collective 477 478 Input Parameter: 479 + ts - timestepping context 480 - arktype - type of ARK-IMEX scheme 481 482 Level: intermediate 483 484 .seealso: TSARKIMEXGetType() 485 @*/ 486 PetscErrorCode TSARKIMEXSetType(TS ts,const TSARKIMEXType arktype) 487 { 488 PetscErrorCode ierr; 489 490 PetscFunctionBegin; 491 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 492 ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,const TSARKIMEXType),(ts,arktype));CHKERRQ(ierr); 493 PetscFunctionReturn(0); 494 } 495 496 #undef __FUNCT__ 497 #define __FUNCT__ "TSARKIMEXGetType" 498 /*@C 499 TSARKIMEXGetType - Get the type of ARK IMEX scheme 500 501 Logically collective 502 503 Input Parameter: 504 . ts - timestepping context 505 506 Output Parameter: 507 . arktype - type of ARK-IMEX scheme 508 509 Level: intermediate 510 511 .seealso: TSARKIMEXGetType() 512 @*/ 513 PetscErrorCode TSARKIMEXGetType(TS ts,const TSARKIMEXType *arktype) 514 { 515 PetscErrorCode ierr; 516 517 PetscFunctionBegin; 518 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 519 ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,const TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr); 520 PetscFunctionReturn(0); 521 } 522 523 EXTERN_C_BEGIN 524 #undef __FUNCT__ 525 #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX" 526 PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,const TSARKIMEXType *arktype) 527 { 528 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 529 PetscErrorCode ierr; 530 531 PetscFunctionBegin; 532 if (!ark->tableau) {ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr);} 533 *arktype = ark->tableau->name; 534 PetscFunctionReturn(0); 535 } 536 #undef __FUNCT__ 537 #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX" 538 PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,const TSARKIMEXType arktype) 539 { 540 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 541 PetscErrorCode ierr; 542 PetscBool match; 543 ARKTableauLink link; 544 545 PetscFunctionBegin; 546 if (ark->tableau) { 547 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 548 if (match) PetscFunctionReturn(0); 549 } 550 for (link = ARKTableauList; link; link=link->next) { 551 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 552 if (match) { 553 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 554 ark->tableau = &link->tab; 555 PetscFunctionReturn(0); 556 } 557 } 558 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 559 PetscFunctionReturn(0); 560 } 561 EXTERN_C_END 562 563 /* ------------------------------------------------------------ */ 564 /*MC 565 TSARKIMEX - ODE solver using Additive Runge-Kutta IMEX schemes 566 567 Level: beginner 568 569 .seealso: TSCreate(), TS, TSSetType() 570 571 M*/ 572 EXTERN_C_BEGIN 573 #undef __FUNCT__ 574 #define __FUNCT__ "TSCreate_ARKIMEX" 575 PetscErrorCode TSCreate_ARKIMEX(TS ts) 576 { 577 TS_ARKIMEX *th; 578 PetscErrorCode ierr; 579 580 PetscFunctionBegin; 581 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 582 ierr = TSARKIMEXInitializePackage(PETSC_NULL);CHKERRQ(ierr); 583 #endif 584 585 ts->ops->reset = TSReset_ARKIMEX; 586 ts->ops->destroy = TSDestroy_ARKIMEX; 587 ts->ops->view = TSView_ARKIMEX; 588 ts->ops->setup = TSSetUp_ARKIMEX; 589 ts->ops->step = TSStep_ARKIMEX; 590 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 591 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 592 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 593 594 ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr); 595 ts->data = (void*)th; 596 597 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","TSARKIMEXGetType_ARKIMEX",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr); 598 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","TSARKIMEXSetType_ARKIMEX",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr); 599 PetscFunctionReturn(0); 600 } 601 EXTERN_C_END 602