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 PetscFunctionReturn(0); 78 } 79 80 #undef __FUNCT__ 81 #define __FUNCT__ "TSARKIMEXRegisterDestroy" 82 /*@C 83 TSARKIMEXRegisterDestroy - Frees the list of schemes that were registered by TSARKIMEXRegister(). 84 85 Not Collective 86 87 Level: advanced 88 89 .keywords: TSARKIMEX, register, destroy 90 .seealso: TSARKIMEXRegister(), TSARKIMEXRegisterAll(), TSARKIMEXRegisterDynamic() 91 @*/ 92 PetscErrorCode TSARKIMEXRegisterDestroy(void) 93 { 94 PetscErrorCode ierr; 95 ARKTableauLink link; 96 97 PetscFunctionBegin; 98 while ((link = ARKTableauList)) { 99 ARKTableau t = &link->tab; 100 ARKTableauList = link->next; 101 ierr = PetscFree6(t->At,t->bt,t->ct,t->A,t->b,t->c);CHKERRQ(ierr); 102 ierr = PetscFree(t->name);CHKERRQ(ierr); 103 ierr = PetscFree(link);CHKERRQ(ierr); 104 } 105 TSARKIMEXRegisterAllCalled = PETSC_FALSE; 106 PetscFunctionReturn(0); 107 } 108 109 #undef __FUNCT__ 110 #define __FUNCT__ "TSARKIMEXInitializePackage" 111 /*@C 112 TSARKIMEXInitializePackage - This function initializes everything in the TSARKIMEX package. It is called 113 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_ARKIMEX() 114 when using static libraries. 115 116 Input Parameter: 117 path - The dynamic library path, or PETSC_NULL 118 119 Level: developer 120 121 .keywords: TS, TSARKIMEX, initialize, package 122 .seealso: PetscInitialize() 123 @*/ 124 PetscErrorCode TSARKIMEXInitializePackage(const char path[]) 125 { 126 PetscErrorCode ierr; 127 128 PetscFunctionBegin; 129 if (TSARKIMEXPackageInitialized) PetscFunctionReturn(0); 130 TSARKIMEXPackageInitialized = PETSC_TRUE; 131 ierr = TSARKIMEXRegisterAll();CHKERRQ(ierr); 132 ierr = PetscRegisterFinalize(TSARKIMEXFinalizePackage);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNCT__ 137 #define __FUNCT__ "TSARKIMEXFinalizePackage" 138 /*@C 139 TSARKIMEXFinalizePackage - This function destroys everything in the TSARKIMEX package. It is 140 called from PetscFinalize(). 141 142 Level: developer 143 144 .keywords: Petsc, destroy, package 145 .seealso: PetscFinalize() 146 @*/ 147 PetscErrorCode TSARKIMEXFinalizePackage(void) 148 { 149 PetscErrorCode ierr; 150 151 PetscFunctionBegin; 152 TSARKIMEXPackageInitialized = PETSC_FALSE; 153 ierr = TSARKIMEXRegisterDestroy();CHKERRQ(ierr); 154 PetscFunctionReturn(0); 155 } 156 157 #undef __FUNCT__ 158 #define __FUNCT__ "TSARKIMEXRegister" 159 PetscErrorCode TSARKIMEXRegister(const TSARKIMEXType name,PetscInt order,PetscInt s, 160 const PetscReal At[],const PetscReal bt[],const PetscReal ct[], 161 const PetscReal A[],const PetscReal b[],const PetscReal c[]) 162 { 163 PetscErrorCode ierr; 164 ARKTableauLink link; 165 ARKTableau t; 166 PetscInt i,j; 167 168 PetscFunctionBegin; 169 ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 170 t = &link->tab; 171 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 172 t->order = order; 173 t->s = s; 174 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); 175 ierr = PetscMemcpy(t->At,At,s*s*sizeof(At[0]));CHKERRQ(ierr); 176 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 177 if (bt) {ierr = PetscMemcpy(t->bt,bt,s*sizeof(bt[0]));CHKERRQ(ierr);} 178 else for (i=0; i<s; i++) t->bt[i] = At[(s-1)*s+i]; 179 if (b) {ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr);} 180 else for (i=0; i<s; i++) t->b[i] = At[(s-1)*s+i]; 181 if (ct) {ierr = PetscMemcpy(t->ct,ct,s*sizeof(ct[0]));CHKERRQ(ierr);} 182 else for (i=0; i<s; i++) for (j=0,t->ct[i]=0; j<s; j++) t->ct[i] += At[i*s+j]; 183 if (c) {ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr);} 184 else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j]; 185 link->next = ARKTableauList; 186 ARKTableauList = link; 187 PetscFunctionReturn(0); 188 } 189 190 #undef __FUNCT__ 191 #define __FUNCT__ "TSStep_ARKIMEX" 192 static PetscErrorCode TSStep_ARKIMEX(TS ts) 193 { 194 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 195 ARKTableau tab = ark->tableau; 196 const PetscInt s = tab->s; 197 const PetscReal *At = tab->At,*A = tab->A,*bt = tab->bt,*b = tab->b,*ct = tab->ct,*c = tab->c; 198 PetscReal *w = ark->work; 199 Vec *Y = ark->Y,*YdotI = ark->YdotI,*YdotRHS = ark->YdotRHS,Ydot = ark->Ydot,W = ark->Work,Z = ark->Z; 200 SNES snes; 201 PetscInt i,j,its,lits; 202 PetscReal h,t; 203 PetscErrorCode ierr; 204 205 PetscFunctionBegin; 206 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 207 h = ts->time_step = ts->next_time_step; 208 t = ts->ptime; 209 210 for (i=0; i<s; i++) { 211 if (At[i*s+i] == 0) { /* This stage is explicit */ 212 ierr = VecCopy(ts->vec_sol,Y[i]);CHKERRQ(ierr); 213 for (j=0; j<i; j++) w[j] = -h*At[i*s+j]; 214 ierr = VecMAXPY(Y[i],i,w,YdotI);CHKERRQ(ierr); 215 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 216 ierr = VecMAXPY(Y[i],i,w,YdotRHS);CHKERRQ(ierr); 217 } else { 218 ark->stage_time = t + h*ct[i]; 219 ark->shift = 1./(h*At[i*s+i]); 220 /* Affine part */ 221 ierr = VecZeroEntries(W);CHKERRQ(ierr); 222 for (j=0; j<i; j++) w[j] = h*A[i*s+j]; 223 ierr = VecMAXPY(W,i,w,YdotRHS);CHKERRQ(ierr); 224 /* Ydot = shift*(Y-Z) */ 225 ierr = VecCopy(ts->vec_sol,Z);CHKERRQ(ierr); 226 for (j=0; j<i; j++) w[j] = h*At[i*s+j]; 227 ierr = VecMAXPY(Z,i,w,YdotRHS);CHKERRQ(ierr); 228 /* Initial guess taken from last stage */ 229 ierr = VecCopy(i>0?Y[i-1]:ts->vec_sol,Y[i]);CHKERRQ(ierr); 230 ierr = SNESSolve(snes,W,Y[i]);CHKERRQ(ierr); 231 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 232 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 233 ts->nonlinear_its += its; ts->linear_its += lits; 234 } 235 ierr = VecZeroEntries(Ydot);CHKERRQ(ierr); 236 ierr = TSComputeIFunction(ts,t+h*ct[i],Y[i],Ydot,YdotI[i],PETSC_TRUE);CHKERRQ(ierr); 237 ierr = TSComputeRHSFunction(ts,t+h*c[i],Y[i],YdotRHS[i]);CHKERRQ(ierr); 238 } 239 for (j=0; j<s; j++) w[j] = h*bt[j]; 240 ierr = VecMAXPY(ts->vec_sol,s,w,YdotI);CHKERRQ(ierr); 241 for (j=0; j<s; j++) w[j] = h*b[j]; 242 ierr = VecMAXPY(ts->vec_sol,s,w,YdotRHS);CHKERRQ(ierr); 243 244 ts->ptime += ts->time_step; 245 ts->next_time_step = ts->time_step; 246 ts->steps++; 247 PetscFunctionReturn(0); 248 } 249 250 /*------------------------------------------------------------*/ 251 #undef __FUNCT__ 252 #define __FUNCT__ "TSReset_ARKIMEX" 253 static PetscErrorCode TSReset_ARKIMEX(TS ts) 254 { 255 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 256 PetscInt s; 257 PetscErrorCode ierr; 258 259 PetscFunctionBegin; 260 if (!ark->tableau) PetscFunctionReturn(0); 261 s = ark->tableau->s; 262 ierr = VecDestroyVecs(s,&ark->Y);CHKERRQ(ierr); 263 ierr = VecDestroyVecs(s,&ark->YdotI);CHKERRQ(ierr); 264 ierr = VecDestroyVecs(s,&ark->YdotRHS);CHKERRQ(ierr); 265 ierr = VecDestroy(&ark->Ydot);CHKERRQ(ierr); 266 ierr = VecDestroy(&ark->Work);CHKERRQ(ierr); 267 ierr = VecDestroy(&ark->Z);CHKERRQ(ierr); 268 ierr = PetscFree(ark->work);CHKERRQ(ierr); 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "TSDestroy_ARKIMEX" 274 static PetscErrorCode TSDestroy_ARKIMEX(TS ts) 275 { 276 PetscErrorCode ierr; 277 278 PetscFunctionBegin; 279 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 280 ierr = PetscFree(ts->data);CHKERRQ(ierr); 281 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","",PETSC_NULL);CHKERRQ(ierr); 282 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","",PETSC_NULL);CHKERRQ(ierr); 283 PetscFunctionReturn(0); 284 } 285 286 /* 287 This defines the nonlinear equation that is to be solved with SNES 288 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 289 */ 290 #undef __FUNCT__ 291 #define __FUNCT__ "SNESTSFormFunction_ARKIMEX" 292 static PetscErrorCode SNESTSFormFunction_ARKIMEX(SNES snes,Vec X,Vec F,TS ts) 293 { 294 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 295 PetscErrorCode ierr; 296 297 PetscFunctionBegin; 298 ierr = VecAXPBYPCZ(ark->Ydot,-ark->shift,ark->shift,0,ark->Z,X);CHKERRQ(ierr); /* Ydot = shift*(X-Z) */ 299 ierr = TSComputeIFunction(ts,ark->stage_time,X,ark->Ydot,X,PETSC_TRUE);CHKERRQ(ierr); 300 PetscFunctionReturn(0); 301 } 302 303 #undef __FUNCT__ 304 #define __FUNCT__ "SNESTSFormJacobian_ARKIMEX" 305 static PetscErrorCode SNESTSFormJacobian_ARKIMEX(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 306 { 307 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 308 PetscErrorCode ierr; 309 310 PetscFunctionBegin; 311 /* ark->Ydot has already been computed in SNESTSFormFunction_ARKIMEX (SNES guarantees this) */ 312 ierr = TSComputeIJacobian(ts,ark->stage_time,X,ark->Ydot,ark->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 #undef __FUNCT__ 317 #define __FUNCT__ "TSSetUp_ARKIMEX" 318 static PetscErrorCode TSSetUp_ARKIMEX(TS ts) 319 { 320 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 321 ARKTableau tab = ark->tableau; 322 PetscInt s = tab->s; 323 PetscErrorCode ierr; 324 325 PetscFunctionBegin; 326 if (!ark->tableau) { 327 ierr = TSSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr); 328 } 329 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->Y);CHKERRQ(ierr); 330 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotI);CHKERRQ(ierr); 331 ierr = VecDuplicateVecs(ts->vec_sol,s,&ark->YdotRHS);CHKERRQ(ierr); 332 ierr = VecDuplicate(ts->vec_sol,&ark->Ydot);CHKERRQ(ierr); 333 ierr = VecDuplicate(ts->vec_sol,&ark->Work);CHKERRQ(ierr); 334 ierr = VecDuplicate(ts->vec_sol,&ark->Z);CHKERRQ(ierr); 335 ierr = PetscMalloc(s*sizeof(ark->work[0]),&ark->work);CHKERRQ(ierr); 336 PetscFunctionReturn(0); 337 } 338 /*------------------------------------------------------------*/ 339 340 #undef __FUNCT__ 341 #define __FUNCT__ "TSSetFromOptions_ARKIMEX" 342 static PetscErrorCode TSSetFromOptions_ARKIMEX(TS ts) 343 { 344 PetscErrorCode ierr; 345 char arktype[256]; 346 347 PetscFunctionBegin; 348 ierr = PetscOptionsHead("ARKIMEX ODE solver options");CHKERRQ(ierr); 349 { 350 ARKTableauLink link; 351 PetscInt count,choice; 352 PetscBool flg; 353 const char **namelist; 354 ierr = PetscStrncpy(arktype,TSARKIMEXDefault,sizeof arktype);CHKERRQ(ierr); 355 for (link=ARKTableauList,count=0; link; link=link->next,count++) ; 356 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 357 for (link=ARKTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 358 ierr = PetscOptionsEList("-ts_arkimex_type","Family of ARK IMEX method","TSARKIMEXSetType",(const char*const*)namelist,count,arktype,&choice,&flg);CHKERRQ(ierr); 359 ierr = TSARKIMEXSetType(ts,flg ? namelist[choice] : arktype);CHKERRQ(ierr); 360 ierr = PetscFree(namelist);CHKERRQ(ierr); 361 } 362 ierr = PetscOptionsTail();CHKERRQ(ierr); 363 PetscFunctionReturn(0); 364 } 365 366 #undef __FUNCT__ 367 #define __FUNCT__ "PetscFormatRealArray" 368 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 369 { 370 int i,left,count; 371 char *p; 372 373 PetscFunctionBegin; 374 for (i=0,p=buf,left=(int)len; i<n; i++) { 375 count = snprintf(p,left,fmt,x[i]); 376 if (count < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"snprintf()"); 377 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 378 left -= count; 379 p += count; 380 *p++ = ' '; 381 } 382 p[i ? 0 : -1] = 0; 383 PetscFunctionReturn(0); 384 } 385 386 #undef __FUNCT__ 387 #define __FUNCT__ "TSView_ARKIMEX" 388 static PetscErrorCode TSView_ARKIMEX(TS ts,PetscViewer viewer) 389 { 390 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 391 ARKTableau tab = ark->tableau; 392 PetscBool iascii; 393 PetscErrorCode ierr; 394 395 PetscFunctionBegin; 396 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 397 if (iascii) { 398 const TSARKIMEXType arktype; 399 char buf[512]; 400 ierr = TSARKIMEXGetType(ts,&arktype);CHKERRQ(ierr); 401 ierr = PetscViewerASCIIPrintf(viewer," ARK IMEX %s\n",arktype);CHKERRQ(ierr); 402 ierr = PetscFormatRealArray(buf,sizeof buf,"%8.6f",tab->s,tab->ct);CHKERRQ(ierr); 403 ierr = PetscViewerASCIIPrintf(viewer," Stiff abscissa at ct = %s\n",buf);CHKERRQ(ierr); 404 } else { 405 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for TS_ARKIMEX",((PetscObject)viewer)->type_name); 406 } 407 PetscFunctionReturn(0); 408 } 409 410 #undef __FUNCT__ 411 #define __FUNCT__ "TSARKIMEXSetType" 412 /*@C 413 TSARKIMEXSetType - Set the type of ARK IMEX scheme 414 415 Logically collective 416 417 Input Parameter: 418 + ts - timestepping context 419 - arktype - type of ARK-IMEX scheme 420 421 Level: intermediate 422 423 .seealso: TSARKIMEXGetType() 424 @*/ 425 PetscErrorCode TSARKIMEXSetType(TS ts,const TSARKIMEXType arktype) 426 { 427 PetscErrorCode ierr; 428 429 PetscFunctionBegin; 430 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 431 ierr = PetscTryMethod(ts,"TSARKIMEXSetType_C",(TS,const TSARKIMEXType),(ts,arktype));CHKERRQ(ierr); 432 PetscFunctionReturn(0); 433 } 434 435 #undef __FUNCT__ 436 #define __FUNCT__ "TSARKIMEXGetType" 437 /*@C 438 TSARKIMEXGetType - Get the type of ARK IMEX scheme 439 440 Logically collective 441 442 Input Parameter: 443 . ts - timestepping context 444 445 Output Parameter: 446 . arktype - type of ARK-IMEX scheme 447 448 Level: intermediate 449 450 .seealso: TSARKIMEXGetType() 451 @*/ 452 PetscErrorCode TSARKIMEXGetType(TS ts,const TSARKIMEXType *arktype) 453 { 454 PetscErrorCode ierr; 455 456 PetscFunctionBegin; 457 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 458 ierr = PetscUseMethod(ts,"TSARKIMEXGetType_C",(TS,const TSARKIMEXType*),(ts,arktype));CHKERRQ(ierr); 459 PetscFunctionReturn(0); 460 } 461 462 EXTERN_C_BEGIN 463 #undef __FUNCT__ 464 #define __FUNCT__ "TSARKIMEXGetType_ARKIMEX" 465 PetscErrorCode TSARKIMEXGetType_ARKIMEX(TS ts,const TSARKIMEXType *arktype) 466 { 467 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 468 PetscErrorCode ierr; 469 470 PetscFunctionBegin; 471 if (!ark->tableau) {ierr = TSARKIMEXSetType(ts,TSARKIMEXDefault);CHKERRQ(ierr);} 472 *arktype = ark->tableau->name; 473 PetscFunctionReturn(0); 474 } 475 #undef __FUNCT__ 476 #define __FUNCT__ "TSARKIMEXSetType_ARKIMEX" 477 PetscErrorCode TSARKIMEXSetType_ARKIMEX(TS ts,const TSARKIMEXType arktype) 478 { 479 TS_ARKIMEX *ark = (TS_ARKIMEX*)ts->data; 480 PetscErrorCode ierr; 481 PetscBool match; 482 ARKTableauLink link; 483 484 PetscFunctionBegin; 485 if (ark->tableau) { 486 ierr = PetscStrcmp(ark->tableau->name,arktype,&match);CHKERRQ(ierr); 487 if (match) PetscFunctionReturn(0); 488 } 489 for (link = ARKTableauList; link; link=link->next) { 490 ierr = PetscStrcmp(link->tab.name,arktype,&match);CHKERRQ(ierr); 491 if (match) { 492 ierr = TSReset_ARKIMEX(ts);CHKERRQ(ierr); 493 ark->tableau = &link->tab; 494 PetscFunctionReturn(0); 495 } 496 } 497 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",arktype); 498 PetscFunctionReturn(0); 499 } 500 EXTERN_C_END 501 502 /* ------------------------------------------------------------ */ 503 /*MC 504 TSARKIMEX - ODE solver using Additive Runge-Kutta IMEX schemes 505 506 Level: beginner 507 508 .seealso: TSCreate(), TS, TSSetType() 509 510 M*/ 511 EXTERN_C_BEGIN 512 #undef __FUNCT__ 513 #define __FUNCT__ "TSCreate_ARKIMEX" 514 PetscErrorCode TSCreate_ARKIMEX(TS ts) 515 { 516 TS_ARKIMEX *th; 517 PetscErrorCode ierr; 518 519 PetscFunctionBegin; 520 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 521 ierr = TSARKIMEXInitializePackage(PETSC_NULL);CHKERRQ(ierr); 522 #endif 523 524 ts->ops->reset = TSReset_ARKIMEX; 525 ts->ops->destroy = TSDestroy_ARKIMEX; 526 ts->ops->view = TSView_ARKIMEX; 527 ts->ops->setup = TSSetUp_ARKIMEX; 528 ts->ops->step = TSStep_ARKIMEX; 529 ts->ops->setfromoptions = TSSetFromOptions_ARKIMEX; 530 ts->ops->snesfunction = SNESTSFormFunction_ARKIMEX; 531 ts->ops->snesjacobian = SNESTSFormJacobian_ARKIMEX; 532 533 ierr = PetscNewLog(ts,TS_ARKIMEX,&th);CHKERRQ(ierr); 534 ts->data = (void*)th; 535 536 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXGetType_C","TSARKIMEXGetType_ARKIMEX",TSARKIMEXGetType_ARKIMEX);CHKERRQ(ierr); 537 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSARKIMEXSetType_C","TSARKIMEXSetType_ARKIMEX",TSARKIMEXSetType_ARKIMEX);CHKERRQ(ierr); 538 PetscFunctionReturn(0); 539 } 540 EXTERN_C_END 541