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 PetscScalar *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 } 466 PetscFunctionReturn(0); 467 } 468 469 #undef __FUNCT__ 470 #define __FUNCT__ "TSARKIMEXSetType" 471 /*@C 472 TSARKIMEXSetType - Set the type of ARK IMEX scheme 473 474 Logically collective 475 476 Input Parameter: 477 + ts - timestepping context 478 - arktype - type of ARK-IMEX scheme 479 480 Level: intermediate 481 482 .seealso: TSARKIMEXGetType() 483 @*/ 484 PetscErrorCode TSARKIMEXSetType(TS ts,const TSARKIMEXType arktype) 485 { 486 PetscErrorCode ierr; 487 488 PetscFunctionBegin; 489 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 490 ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,const TSARKIMEXType),(ts,arktype));CHKERRQ(ierr); 491 PetscFunctionReturn(0); 492 } 493 494 #undef __FUNCT__ 495 #define __FUNCT__ "TSARKIMEXGetType" 496 /*@C 497 TSARKIMEXGetType - Get the type of ARK IMEX scheme 498 499 Logically collective 500 501 Input Parameter: 502 . ts - timestepping context 503 504 Output Parameter: 505 . arktype - type of ARK-IMEX scheme 506 507 Level: intermediate 508 509 .seealso: TSARKIMEXGetType() 510 @*/ 511 PetscErrorCode TSARKIMEXGetType(TS ts,const TSARKIMEXType *arktype) 512 { 513 PetscErrorCode ierr; 514 515 PetscFunctionBegin; 516 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 517 ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,const TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr); 518 PetscFunctionReturn(0); 519 } 520 521 EXTERN_C_BEGIN 522 #undef __FUNCT__ 523 #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX" 524 PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,const TSARKIMEXType *arktype) 525 { 526 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 527 PetscErrorCode ierr; 528 529 PetscFunctionBegin; 530 if (!ark->tableau) {ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr);} 531 *arktype = ark->tableau->name; 532 PetscFunctionReturn(0); 533 } 534 #undef __FUNCT__ 535 #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX" 536 PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,const TSARKIMEXType arktype) 537 { 538 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 539 PetscErrorCode ierr; 540 PetscBool match; 541 ARKTableauLink link; 542 543 PetscFunctionBegin; 544 if (ark->tableau) { 545 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 546 if (match) PetscFunctionReturn(0); 547 } 548 for (link = ARKTableauList; link; link=link->next) { 549 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 550 if (match) { 551 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 552 ark->tableau = &link->tab; 553 PetscFunctionReturn(0); 554 } 555 } 556 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 557 PetscFunctionReturn(0); 558 } 559 EXTERN_C_END 560 561 /* ------------------------------------------------------------ */ 562 /*MC 563 TSARKIMEX - ODE solver using Additive Runge-Kutta IMEX schemes 564 565 Level: beginner 566 567 .seealso: TSCreate(), TS, TSSetType() 568 569 M*/ 570 EXTERN_C_BEGIN 571 #undef __FUNCT__ 572 #define __FUNCT__ "TSCreate_ARKIMEX" 573 PetscErrorCode TSCreate_ARKIMEX(TS ts) 574 { 575 TS_ARKIMEX *th; 576 PetscErrorCode ierr; 577 578 PetscFunctionBegin; 579 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 580 ierr = TSARKIMEXInitializePackage(PETSC_NULL);CHKERRQ(ierr); 581 #endif 582 583 ts->ops->reset = TSReset_ARKIMEX; 584 ts->ops->destroy = TSDestroy_ARKIMEX; 585 ts->ops->view = TSView_ARKIMEX; 586 ts->ops->setup = TSSetUp_ARKIMEX; 587 ts->ops->step = TSStep_ARKIMEX; 588 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 589 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 590 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 591 592 ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr); 593 ts->data = (void*)th; 594 595 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","TSARKIMEXGetType_ARKIMEX",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr); 596 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","TSARKIMEXSetType_ARKIMEX",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr); 597 PetscFunctionReturn(0); 598 } 599 EXTERN_C_END 600