1 /* 2 Code for timestepping with Rosenbrock W methods 3 4 Notes: 5 The general system is written as 6 7 G(t,X,Xdot) = F(t,X) 8 9 where G represents the stiff part of the physics and F represents the non-stiff part. 10 This method is designed to be linearly implicit on G and can use an approximate and lagged Jacobian. 11 12 */ 13 #include <private/tsimpl.h> /*I "petscts.h" I*/ 14 15 #include <../src/mat/blockinvert.h> 16 17 static const TSRosWType TSRosWDefault = TSROSW2P; 18 static PetscBool TSRosWRegisterAllCalled; 19 static PetscBool TSRosWPackageInitialized; 20 21 typedef struct _RosWTableau *RosWTableau; 22 struct _RosWTableau { 23 char *name; 24 PetscInt order; /* Classical approximation order of the method */ 25 PetscInt s; /* Number of stages */ 26 PetscReal *A; /* Propagation table, strictly lower triangular */ 27 PetscReal *Gamma; /* Stage table, lower triangular with nonzero diagonal */ 28 PetscBool *GammaZeroDiag; /* Diagonal entries that are zero in stage table Gamma, vector indicating explicit statages */ 29 PetscReal *b; /* Step completion table */ 30 PetscReal *bembed; /* Step completion table for embedded method of order one less */ 31 PetscReal *ASum; /* Row sum of A */ 32 PetscReal *GammaSum; /* Row sum of Gamma, only needed for non-autonomous systems */ 33 PetscReal *At; /* Propagation table in transformed variables */ 34 PetscReal *bt; /* Step completion table in transformed variables */ 35 PetscReal *bembedt; /* Step completion table of order one less in transformed variables */ 36 PetscReal *GammaInv; /* Inverse of Gamma, used for transformed variables */ 37 PetscReal ccfl; /* Placeholder for CFL coefficient relative to forward Euler */ 38 }; 39 typedef struct _RosWTableauLink *RosWTableauLink; 40 struct _RosWTableauLink { 41 struct _RosWTableau tab; 42 RosWTableauLink next; 43 }; 44 static RosWTableauLink RosWTableauList; 45 46 typedef struct { 47 RosWTableau tableau; 48 Vec *Y; /* States computed during the step, used to complete the step */ 49 Vec Ydot; /* Work vector holding Ydot during residual evaluation */ 50 Vec Ystage; /* Work vector for the state value at each stage */ 51 Vec Zdot; /* Ydot = Zdot + shift*Y */ 52 Vec Zstage; /* Y = Zstage + Y */ 53 PetscScalar *work; /* Scalar work space of length number of stages, used to prepare VecMAXPY() */ 54 PetscReal shift; 55 PetscReal stage_time; 56 PetscReal stage_explicit; /* Flag indicates that the current stage is explicit */ 57 PetscBool recompute_jacobian; /* Recompute the Jacobian at each stage, default is to freeze the Jacobian at the start of each step */ 58 PetscBool step_taken; /* ts->vec_sol has been advanced to the end of the current time step */ 59 } TS_RosW; 60 61 /*MC 62 TSROSW2M - Two stage second order L-stable Rosenbrock-W scheme. 63 64 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2P. 65 66 Level: intermediate 67 68 .seealso: TSROSW 69 M*/ 70 71 /*MC 72 TSROSW2P - Two stage second order L-stable Rosenbrock-W scheme. 73 74 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. This method is a reflection of TSROSW2M. 75 76 Level: intermediate 77 78 .seealso: TSROSW 79 M*/ 80 81 /*MC 82 TSROSWRA3PW - Three stage third order Rosenbrock-W scheme for PDAE of index 1. 83 84 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 85 86 This is strongly A-stable with R(infty) = 0.73. The embedded method of order 2 is strongly A-stable with R(infty) = 0.73. 87 88 References: 89 Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005. 90 91 Level: intermediate 92 93 .seealso: TSROSW 94 M*/ 95 96 /*MC 97 TSROSWRA34PW2 - Four stage third order L-stable Rosenbrock-W scheme for PDAE of index 1. 98 99 Only an approximate Jacobian is needed. By default, it is only recomputed once per step. 100 101 This is strongly A-stable with R(infty) = 0. The embedded method of order 2 is strongly A-stable with R(infty) = 0.48. 102 103 References: 104 Rang and Angermann, New Rosenbrock-W methods of order 3 for partial differential algebraic equations of index 1, 2005. 105 106 Level: intermediate 107 108 .seealso: TSROSW 109 M*/ 110 111 /*MC 112 TSROSWRODAS3 - Four stage third order L-stable Rosenbrock scheme 113 114 By default, the Jacobian is only recomputed once per step. 115 116 Both the third order and embedded second order methods are stiffly accurate and L-stable. 117 118 References: 119 Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 120 121 Level: intermediate 122 123 .seealso: TSROSW, TSROSWSANDU3 124 M*/ 125 126 /*MC 127 TSROSWSANDU3 - Three stage third order L-stable Rosenbrock scheme 128 129 By default, the Jacobian is only recomputed once per step. 130 131 The third order method is L-stable, but not stiffly accurate. 132 The second order embedded method is strongly A-stable with R(infty) = 0.5. 133 The internal stages are L-stable. 134 This method is called ROS3 in the paper. 135 136 References: 137 Sandu et al, Benchmarking stiff ODE solvers for atmospheric chemistry problems II, Rosenbrock solvers, 1997. 138 139 Level: intermediate 140 141 .seealso: TSROSW, TSROSWRODAS3 142 M*/ 143 144 /*MC 145 TSROSWASSP3P3S1C - A-stable Rosenbrock-W method with SSP explicit part, third order, three stages 146 147 By default, the Jacobian is only recomputed once per step. 148 149 A-stable SPP explicit order 3, 3 stages, CFL 1 (eff = 1/3) 150 151 References: 152 Emil Constantinescu 153 154 Level: intermediate 155 156 .seealso: TSROSW, TSROSWLASSP3P4S2C, TSROSWLLSSP3P3S2C, SSP 157 M*/ 158 159 /*MC 160 TSROSWLASSP3P4S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, three stages 161 162 By default, the Jacobian is only recomputed once per step. 163 164 L-stable (A-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 165 166 References: 167 Emil Constantinescu 168 169 Level: intermediate 170 171 .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLLSSP3P3S2C, TSSSP 172 M*/ 173 174 /*MC 175 TSROSWLLSSP3P3S2C - L-stable Rosenbrock-W method with SSP explicit part, third order, three stages 176 177 By default, the Jacobian is only recomputed once per step. 178 179 L-stable (L-stable embedded) SPP explicit order 3, 4 stages, CFL 2 (eff = 1/2) 180 181 References: 182 Emil Constantinescu 183 184 Level: intermediate 185 186 .seealso: TSROSW, TSROSWASSP3P3S1C, TSROSWLASSP3P4S2C, TSSSP 187 M*/ 188 189 #undef __FUNCT__ 190 #define __FUNCT__ "TSRosWRegisterAll" 191 /*@C 192 TSRosWRegisterAll - Registers all of the additive Runge-Kutta implicit-explicit methods in TSRosW 193 194 Not Collective, but should be called by all processes which will need the schemes to be registered 195 196 Level: advanced 197 198 .keywords: TS, TSRosW, register, all 199 200 .seealso: TSRosWRegisterDestroy() 201 @*/ 202 PetscErrorCode TSRosWRegisterAll(void) 203 { 204 PetscErrorCode ierr; 205 206 PetscFunctionBegin; 207 if (TSRosWRegisterAllCalled) PetscFunctionReturn(0); 208 TSRosWRegisterAllCalled = PETSC_TRUE; 209 210 { 211 const PetscReal g = 1. + 1./PetscSqrtReal(2.0); 212 const PetscReal 213 A[2][2] = {{0,0}, {1.,0}}, 214 Gamma[2][2] = {{g,0}, {-2.*g,g}}, 215 b[2] = {0.5,0.5}, 216 b1[2] = {1.0,0.0}; 217 ierr = TSRosWRegister(TSROSW2P,2,2,&A[0][0],&Gamma[0][0],b,b1);CHKERRQ(ierr); 218 } 219 { 220 const PetscReal g = 1. - 1./PetscSqrtReal(2.0); 221 const PetscReal 222 A[2][2] = {{0,0}, {1.,0}}, 223 Gamma[2][2] = {{g,0}, {-2.*g,g}}, 224 b[2] = {0.5,0.5}, 225 b1[2] = {1.0,0.0}; 226 ierr = TSRosWRegister(TSROSW2M,2,2,&A[0][0],&Gamma[0][0],b,b1);CHKERRQ(ierr); 227 } 228 { 229 const PetscReal g = 7.8867513459481287e-01; 230 const PetscReal 231 A[3][3] = {{0,0,0}, 232 {1.5773502691896257e+00,0,0}, 233 {0.5,0,0}}, 234 Gamma[3][3] = {{g,0,0}, 235 {-1.5773502691896257e+00,g,0}, 236 {-6.7075317547305480e-01,1.7075317547305482e-01,g}}, 237 b[3] = {1.0566243270259355e-01,4.9038105676657971e-02,8.4529946162074843e-01}, 238 b2[3] = {-1.7863279495408180e-01,1./3.,8.4529946162074843e-01}; 239 ierr = TSRosWRegister(TSROSWRA3PW,3,3,&A[0][0],&Gamma[0][0],b,b2);CHKERRQ(ierr); 240 } 241 { 242 const PetscReal g = 4.3586652150845900e-01; 243 const PetscReal 244 A[4][4] = {{0,0,0,0}, 245 {8.7173304301691801e-01,0,0,0}, 246 {8.4457060015369423e-01,-1.1299064236484185e-01,0,0}, 247 {0,0,1.,0}}, 248 Gamma[4][4] = {{g,0,0,0}, 249 {-8.7173304301691801e-01,g,0,0}, 250 {-9.0338057013044082e-01,5.4180672388095326e-02,g,0}, 251 {2.4212380706095346e-01,-1.2232505839045147e+00,5.4526025533510214e-01,g}}, 252 b[4] = {2.4212380706095346e-01,-1.2232505839045147e+00,1.5452602553351020e+00,4.3586652150845900e-01}, 253 b2[4] = {3.7810903145819369e-01,-9.6042292212423178e-02,5.0000000000000000e-01,2.1793326075422950e-01}; 254 ierr = TSRosWRegister(TSROSWRA34PW2,3,4,&A[0][0],&Gamma[0][0],b,b2);CHKERRQ(ierr); 255 } 256 { 257 const PetscReal g = 0.5; 258 const PetscReal 259 A[4][4] = {{0,0,0,0}, 260 {0,0,0,0}, 261 {1.,0,0,0}, 262 {0.75,-0.25,0.5,0}}, 263 Gamma[4][4] = {{g,0,0,0}, 264 {1.,g,0,0}, 265 {-0.25,-0.25,g,0}, 266 {1./12,1./12,-2./3,g}}, 267 b[4] = {5./6,-1./6,-1./6,0.5}, 268 b2[4] = {0.75,-0.25,0.5,0}; 269 ierr = TSRosWRegister(TSROSWRODAS3,3,4,&A[0][0],&Gamma[0][0],b,b2);CHKERRQ(ierr); 270 } 271 { 272 const PetscReal g = 0.43586652150845899941601945119356; 273 const PetscReal 274 A[3][3] = {{0,0,0}, 275 {g,0,0}, 276 {g,0,0}}, 277 Gamma[3][3] = {{g,0,0}, 278 {-0.19294655696029095575009695436041,g,0}, 279 {0,1.74927148125794685173529749738960,g}}, 280 b[3] = {-0.75457412385404315829818998646589,1.94100407061964420292840123379419,-0.18642994676560104463021124732829}, 281 b2[3] = {-1.53358745784149585370766523913002,2.81745131148625772213931745457622,-0.28386385364476186843165221544619}; 282 ierr = TSRosWRegister(TSROSWSANDU3,3,3,&A[0][0],&Gamma[0][0],b,b2);CHKERRQ(ierr); 283 } 284 { 285 const PetscReal g = (3.0+sqrt(3.0))/6.0; 286 const PetscReal 287 A[3][3] = {{0,0,0}, 288 {1,0,0}, 289 {0.25,0.25,0}}, 290 Gamma[3][3] = {{0,0,0}, 291 {(-3.0-sqrt(3.0))/6.0,g,0}, 292 {(-3.0-sqrt(3.0))/24.0,(-3.0-sqrt(3.0))/8.0,g}}, 293 b[3] = {1./6.,1./6.,2./3.}, 294 b2[3] = {1./4.,1./4.,1./2.}; 295 ierr = TSRosWRegister(TSROSWASSP3P3S1C,3,3,&A[0][0],&Gamma[0][0],b,b2);CHKERRQ(ierr); 296 } 297 298 { 299 const PetscReal 300 A[4][4] = {{0,0,0,0}, 301 {1./2.,0,0,0}, 302 {1./2.,1./2.,0,0}, 303 {1./6.,1./6.,1./6.,0}}, 304 Gamma[4][4] = {{1./2.,0,0,0}, 305 {0.0,1./4.,0,0}, 306 {-2.,-2./3.,2./3.,0}, 307 {1./2.,5./36.,-2./9,0}}, 308 b[4] = {1./6.,1./6.,1./6.,1./2.}, 309 b2[4] = {1./8.,3./4.,1./8.,0}; 310 ierr = TSRosWRegister(TSROSWLASSP3P4S2C,3,4,&A[0][0],&Gamma[0][0],b,b2);CHKERRQ(ierr); 311 } 312 313 { 314 const PetscReal 315 A[4][4] = {{0,0,0,0}, 316 {1./2.,0,0,0}, 317 {1./2.,1./2.,0,0}, 318 {1./6.,1./6.,1./6.,0}}, 319 Gamma[4][4] = {{1./2.,0,0,0}, 320 {0.0,3./4.,0,0}, 321 {-2./3.,-23./9.,2./9.,0}, 322 {1./18.,65./108.,-2./27,0}}, 323 b[4] = {1./6.,1./6.,1./6.,1./2.}, 324 b2[4] = {3./16.,10./16.,3./16.,0}; 325 ierr = TSRosWRegister(TSROSWLLSSP3P3S2C,3,4,&A[0][0],&Gamma[0][0],b,b2);CHKERRQ(ierr); 326 } 327 PetscFunctionReturn(0); 328 } 329 330 #undef __FUNCT__ 331 #define __FUNCT__ "TSRosWRegisterDestroy" 332 /*@C 333 TSRosWRegisterDestroy - Frees the list of schemes that were registered by TSRosWRegister(). 334 335 Not Collective 336 337 Level: advanced 338 339 .keywords: TSRosW, register, destroy 340 .seealso: TSRosWRegister(), TSRosWRegisterAll(), TSRosWRegisterDynamic() 341 @*/ 342 PetscErrorCode TSRosWRegisterDestroy(void) 343 { 344 PetscErrorCode ierr; 345 RosWTableauLink link; 346 347 PetscFunctionBegin; 348 while ((link = RosWTableauList)) { 349 RosWTableau t = &link->tab; 350 RosWTableauList = link->next; 351 ierr = PetscFree5(t->A,t->Gamma,t->b,t->ASum,t->GammaSum);CHKERRQ(ierr); 352 ierr = PetscFree4(t->At,t->bt,t->GammaInv,t->GammaZeroDiag);CHKERRQ(ierr); 353 ierr = PetscFree2(t->bembed,t->bembedt);CHKERRQ(ierr); 354 ierr = PetscFree(t->name);CHKERRQ(ierr); 355 ierr = PetscFree(link);CHKERRQ(ierr); 356 } 357 TSRosWRegisterAllCalled = PETSC_FALSE; 358 PetscFunctionReturn(0); 359 } 360 361 #undef __FUNCT__ 362 #define __FUNCT__ "TSRosWInitializePackage" 363 /*@C 364 TSRosWInitializePackage - This function initializes everything in the TSRosW package. It is called 365 from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSCreate_RosW() 366 when using static libraries. 367 368 Input Parameter: 369 path - The dynamic library path, or PETSC_NULL 370 371 Level: developer 372 373 .keywords: TS, TSRosW, initialize, package 374 .seealso: PetscInitialize() 375 @*/ 376 PetscErrorCode TSRosWInitializePackage(const char path[]) 377 { 378 PetscErrorCode ierr; 379 380 PetscFunctionBegin; 381 if (TSRosWPackageInitialized) PetscFunctionReturn(0); 382 TSRosWPackageInitialized = PETSC_TRUE; 383 ierr = TSRosWRegisterAll();CHKERRQ(ierr); 384 ierr = PetscRegisterFinalize(TSRosWFinalizePackage);CHKERRQ(ierr); 385 PetscFunctionReturn(0); 386 } 387 388 #undef __FUNCT__ 389 #define __FUNCT__ "TSRosWFinalizePackage" 390 /*@C 391 TSRosWFinalizePackage - This function destroys everything in the TSRosW package. It is 392 called from PetscFinalize(). 393 394 Level: developer 395 396 .keywords: Petsc, destroy, package 397 .seealso: PetscFinalize() 398 @*/ 399 PetscErrorCode TSRosWFinalizePackage(void) 400 { 401 PetscErrorCode ierr; 402 403 PetscFunctionBegin; 404 TSRosWPackageInitialized = PETSC_FALSE; 405 ierr = TSRosWRegisterDestroy();CHKERRQ(ierr); 406 PetscFunctionReturn(0); 407 } 408 409 #undef __FUNCT__ 410 #define __FUNCT__ "TSRosWRegister" 411 /*@C 412 TSRosWRegister - register a Rosenbrock W scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation 413 414 Not Collective, but the same schemes should be registered on all processes on which they will be used 415 416 Input Parameters: 417 + name - identifier for method 418 . order - approximation order of method 419 . s - number of stages, this is the dimension of the matrices below 420 . A - Table of propagated stage coefficients (dimension s*s, row-major), strictly lower triangular 421 . Gamma - Table of coefficients in implicit stage equations (dimension s*s, row-major), lower triangular with nonzero diagonal 422 . b - Step completion table (dimension s) 423 - bembed - Step completion table for a scheme of order one less (dimension s, PETSC_NULL if no embedded scheme is available) 424 425 Notes: 426 Several Rosenbrock W methods are provided, this function is only needed to create new methods. 427 428 Level: advanced 429 430 .keywords: TS, register 431 432 .seealso: TSRosW 433 @*/ 434 PetscErrorCode TSRosWRegister(const TSRosWType name,PetscInt order,PetscInt s, 435 const PetscReal A[],const PetscReal Gamma[],const PetscReal b[],const PetscReal bembed[]) 436 { 437 PetscErrorCode ierr; 438 RosWTableauLink link; 439 RosWTableau t; 440 PetscInt i,j,k; 441 PetscScalar *GammaInv; 442 443 PetscFunctionBegin; 444 PetscValidCharPointer(name,1); 445 PetscValidPointer(A,4); 446 PetscValidPointer(Gamma,5); 447 PetscValidPointer(b,6); 448 if (bembed) PetscValidPointer(bembed,7); 449 450 ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr); 451 ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr); 452 t = &link->tab; 453 ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr); 454 t->order = order; 455 t->s = s; 456 ierr = PetscMalloc5(s*s,PetscReal,&t->A,s*s,PetscReal,&t->Gamma,s,PetscReal,&t->b,s,PetscReal,&t->ASum,s,PetscReal,&t->GammaSum);CHKERRQ(ierr); 457 ierr = PetscMalloc4(s*s,PetscReal,&t->At,s,PetscReal,&t->bt,s*s,PetscReal,&t->GammaInv,s,PetscBool,&t->GammaZeroDiag);CHKERRQ(ierr); 458 ierr = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr); 459 ierr = PetscMemcpy(t->Gamma,Gamma,s*s*sizeof(Gamma[0]));CHKERRQ(ierr); 460 ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); 461 if (bembed) { 462 ierr = PetscMalloc2(s,PetscReal,&t->bembed,s,PetscReal,&t->bembedt);CHKERRQ(ierr); 463 ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr); 464 } 465 for (i=0; i<s; i++) { 466 t->ASum[i] = 0; 467 t->GammaSum[i] = 0; 468 for (j=0; j<s; j++) { 469 t->ASum[i] += A[i*s+j]; 470 t->GammaSum[i] += Gamma[i*s+j]; 471 } 472 } 473 ierr = PetscMalloc(s*s*sizeof(PetscScalar),&GammaInv);CHKERRQ(ierr); /* Need to use Scalar for inverse, then convert back to Real */ 474 for (i=0; i<s*s; i++) GammaInv[i] = Gamma[i]; 475 for (i=0; i<s; i++) { 476 if (Gamma[i*s+i] == 0.0) { 477 GammaInv[i*s+i] = 1.0; 478 t->GammaZeroDiag[i] = PETSC_TRUE; 479 } else { 480 t->GammaZeroDiag[i] = PETSC_FALSE; 481 } 482 } 483 484 switch (s) { 485 case 1: GammaInv[0] = 1./GammaInv[0]; break; 486 case 2: ierr = Kernel_A_gets_inverse_A_2(GammaInv,0);CHKERRQ(ierr); break; 487 case 3: ierr = Kernel_A_gets_inverse_A_3(GammaInv,0);CHKERRQ(ierr); break; 488 case 4: ierr = Kernel_A_gets_inverse_A_4(GammaInv,0);CHKERRQ(ierr); break; 489 case 5: { 490 PetscInt ipvt5[5]; 491 MatScalar work5[5*5]; 492 ierr = Kernel_A_gets_inverse_A_5(GammaInv,ipvt5,work5,0);CHKERRQ(ierr); break; 493 } 494 case 6: ierr = Kernel_A_gets_inverse_A_6(GammaInv,0);CHKERRQ(ierr); break; 495 case 7: ierr = Kernel_A_gets_inverse_A_7(GammaInv,0);CHKERRQ(ierr); break; 496 default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented for %D stages",s); 497 } 498 for (i=0; i<s*s; i++) t->GammaInv[i] = PetscRealPart(GammaInv[i]); 499 ierr = PetscFree(GammaInv);CHKERRQ(ierr); 500 for (i=0; i<s; i++) { 501 for (j=0; j<s; j++) { 502 t->At[i*s+j] = 0; 503 for (k=0; k<s; k++) { 504 t->At[i*s+j] += t->A[i*s+k] * t->GammaInv[k*s+j]; 505 } 506 } 507 t->bt[i] = 0; 508 for (j=0; j<s; j++) { 509 t->bt[i] += t->b[j] * t->GammaInv[j*s+i]; 510 } 511 if (bembed) { 512 t->bembedt[i] = 0; 513 for (j=0; j<s; j++) { 514 t->bembedt[i] += t->bembed[j] * t->GammaInv[j*s+i]; 515 } 516 } 517 } 518 t->ccfl = 1.0; /* Fix this */ 519 520 link->next = RosWTableauList; 521 RosWTableauList = link; 522 PetscFunctionReturn(0); 523 } 524 525 #undef __FUNCT__ 526 #define __FUNCT__ "TSEvaluateStep_RosW" 527 /* 528 The step completion formula is 529 530 x1 = x0 + b^T Y 531 532 where Y is the multi-vector of stages corrections. This function can be called before or after ts->vec_sol has been 533 updated. Suppose we have a completion formula b and an embedded formula be of different order. We can write 534 535 x1e = x0 + be^T Y 536 = x1 - b^T Y + be^T Y 537 = x1 + (be - b)^T Y 538 539 so we can evaluate the method of different order even after the step has been optimistically completed. 540 */ 541 static PetscErrorCode TSEvaluateStep_RosW(TS ts,PetscInt order,Vec X,PetscBool *done) 542 { 543 TS_RosW *ros = (TS_RosW*)ts->data; 544 RosWTableau tab = ros->tableau; 545 PetscScalar *w = ros->work; 546 PetscInt i; 547 PetscErrorCode ierr; 548 549 PetscFunctionBegin; 550 if (order == tab->order) { 551 if (ros->step_taken) {ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr);} 552 else { 553 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 554 ierr = VecMAXPY(X,tab->s,tab->bt,ros->Y);CHKERRQ(ierr); 555 } 556 if (done) *done = PETSC_TRUE; 557 PetscFunctionReturn(0); 558 } else if (order == tab->order-1) { 559 if (!tab->bembedt) goto unavailable; 560 if (ros->step_taken) { 561 for (i=0; i<tab->s; i++) w[i] = tab->bembedt[i] - tab->bt[i]; 562 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 563 ierr = VecMAXPY(X,tab->s,w,ros->Y);CHKERRQ(ierr); 564 } else { 565 ierr = VecCopy(ts->vec_sol,X);CHKERRQ(ierr); 566 ierr = VecMAXPY(X,tab->s,tab->bembedt,ros->Y);CHKERRQ(ierr); 567 } 568 if (done) *done = PETSC_TRUE; 569 PetscFunctionReturn(0); 570 } 571 unavailable: 572 if (done) *done = PETSC_FALSE; 573 else SETERRQ3(((PetscObject)ts)->comm,PETSC_ERR_SUP,"Rosenbrock-W '%s' of order %D cannot evaluate step at order %D",tab->name,tab->order,order); 574 PetscFunctionReturn(0); 575 } 576 577 #undef __FUNCT__ 578 #define __FUNCT__ "TSStep_RosW" 579 static PetscErrorCode TSStep_RosW(TS ts) 580 { 581 TS_RosW *ros = (TS_RosW*)ts->data; 582 RosWTableau tab = ros->tableau; 583 const PetscInt s = tab->s; 584 const PetscReal *At = tab->At,*Gamma = tab->Gamma,*ASum = tab->ASum,*GammaInv = tab->GammaInv; 585 const PetscBool *GammaZeroDiag = tab->GammaZeroDiag; 586 PetscScalar *w = ros->work; 587 Vec *Y = ros->Y,Zdot = ros->Zdot,Zstage = ros->Zstage; 588 SNES snes; 589 TSAdapt adapt; 590 PetscInt i,j,its,lits,reject,next_scheme; 591 PetscReal next_time_step; 592 PetscBool accept; 593 PetscErrorCode ierr; 594 595 PetscFunctionBegin; 596 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 597 next_time_step = ts->time_step; 598 accept = PETSC_TRUE; 599 ros->step_taken = PETSC_FALSE; 600 601 for (reject=0; reject<ts->max_reject; reject++,ts->reject++) { 602 const PetscReal h = ts->time_step; 603 for (i=0; i<s; i++) { 604 ros->stage_time = ts->ptime + h*ASum[i]; 605 if (GammaZeroDiag[i]) { 606 ros->stage_explicit = PETSC_TRUE; 607 ros->shift = 1./h; 608 } else { 609 ros->stage_explicit = PETSC_FALSE; 610 ros->shift = 1./(h*Gamma[i*s+i]); 611 } 612 613 ierr = VecCopy(ts->vec_sol,Zstage);CHKERRQ(ierr); 614 ierr = VecMAXPY(Zstage,i,&At[i*s+0],Y);CHKERRQ(ierr); 615 616 for (j=0; j<i; j++) w[j] = 1./h * GammaInv[i*s+j]; 617 ierr = VecZeroEntries(Zdot);CHKERRQ(ierr); 618 ierr = VecMAXPY(Zdot,i,w,Y);CHKERRQ(ierr); 619 620 /* Initial guess taken from last stage */ 621 ierr = VecZeroEntries(Y[i]);CHKERRQ(ierr); 622 623 if (!ros->recompute_jacobian && !i) { 624 ierr = SNESSetLagJacobian(snes,-2);CHKERRQ(ierr); /* Recompute the Jacobian on this solve, but not again */ 625 } 626 627 ierr = SNESSolve(snes,PETSC_NULL,Y[i]);CHKERRQ(ierr); 628 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 629 ierr = SNESGetLinearSolveIterations(snes,&lits);CHKERRQ(ierr); 630 ts->nonlinear_its += its; ts->linear_its += lits; 631 } 632 ierr = TSEvaluateStep(ts,tab->order,ts->vec_sol,PETSC_NULL);CHKERRQ(ierr); 633 ros->step_taken = PETSC_TRUE; 634 635 /* Register only the current method as a candidate because we're not supporting multiple candidates yet. */ 636 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 637 ierr = TSAdaptCandidatesClear(adapt);CHKERRQ(ierr); 638 ierr = TSAdaptCandidateAdd(adapt,tab->name,tab->order,1,tab->ccfl,1.*tab->s,PETSC_TRUE);CHKERRQ(ierr); 639 ierr = TSAdaptChoose(adapt,ts,ts->time_step,&next_scheme,&next_time_step,&accept);CHKERRQ(ierr); 640 if (accept) { 641 /* ignore next_scheme for now */ 642 ts->ptime += ts->time_step; 643 ts->time_step = next_time_step; 644 ts->steps++; 645 break; 646 } else { /* Roll back the current step */ 647 for (i=0; i<s; i++) w[i] = -tab->bt[i]; 648 ierr = VecMAXPY(ts->vec_sol,s,w,Y);CHKERRQ(ierr); 649 ts->time_step = next_time_step; 650 ros->step_taken = PETSC_FALSE; 651 } 652 } 653 654 PetscFunctionReturn(0); 655 } 656 657 #undef __FUNCT__ 658 #define __FUNCT__ "TSInterpolate_RosW" 659 static PetscErrorCode TSInterpolate_RosW(TS ts,PetscReal itime,Vec X) 660 { 661 TS_RosW *ros = (TS_RosW*)ts->data; 662 663 PetscFunctionBegin; 664 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_SUP,"TSRosW %s does not have an interpolation formula",ros->tableau->name); 665 PetscFunctionReturn(0); 666 } 667 668 /*------------------------------------------------------------*/ 669 #undef __FUNCT__ 670 #define __FUNCT__ "TSReset_RosW" 671 static PetscErrorCode TSReset_RosW(TS ts) 672 { 673 TS_RosW *ros = (TS_RosW*)ts->data; 674 PetscInt s; 675 PetscErrorCode ierr; 676 677 PetscFunctionBegin; 678 if (!ros->tableau) PetscFunctionReturn(0); 679 s = ros->tableau->s; 680 ierr = VecDestroyVecs(s,&ros->Y);CHKERRQ(ierr); 681 ierr = VecDestroy(&ros->Ydot);CHKERRQ(ierr); 682 ierr = VecDestroy(&ros->Ystage);CHKERRQ(ierr); 683 ierr = VecDestroy(&ros->Zdot);CHKERRQ(ierr); 684 ierr = VecDestroy(&ros->Zstage);CHKERRQ(ierr); 685 ierr = PetscFree(ros->work);CHKERRQ(ierr); 686 PetscFunctionReturn(0); 687 } 688 689 #undef __FUNCT__ 690 #define __FUNCT__ "TSDestroy_RosW" 691 static PetscErrorCode TSDestroy_RosW(TS ts) 692 { 693 PetscErrorCode ierr; 694 695 PetscFunctionBegin; 696 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 697 ierr = PetscFree(ts->data);CHKERRQ(ierr); 698 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","",PETSC_NULL);CHKERRQ(ierr); 699 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","",PETSC_NULL);CHKERRQ(ierr); 700 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","",PETSC_NULL);CHKERRQ(ierr); 701 PetscFunctionReturn(0); 702 } 703 704 /* 705 This defines the nonlinear equation that is to be solved with SNES 706 G(U) = F[t0+Theta*dt, U, (U-U0)*shift] = 0 707 */ 708 #undef __FUNCT__ 709 #define __FUNCT__ "SNESTSFormFunction_RosW" 710 static PetscErrorCode SNESTSFormFunction_RosW(SNES snes,Vec X,Vec F,TS ts) 711 { 712 TS_RosW *ros = (TS_RosW*)ts->data; 713 PetscErrorCode ierr; 714 715 PetscFunctionBegin; 716 if (ros->stage_explicit) { 717 ierr = VecAXPBY(ros->Ydot,ros->shift,0.0,X);CHKERRQ(ierr); /* Ydot = shift*X*/ 718 } else { 719 ierr = VecWAXPY(ros->Ydot,ros->shift,X,ros->Zdot);CHKERRQ(ierr); /* Ydot = shift*X + Zdot */ 720 } 721 ierr = VecWAXPY(ros->Ystage,1.0,X,ros->Zstage);CHKERRQ(ierr); /* Ystage = X + Zstage */ 722 ierr = TSComputeIFunction(ts,ros->stage_time,ros->Ystage,ros->Ydot,F,PETSC_FALSE);CHKERRQ(ierr); 723 PetscFunctionReturn(0); 724 } 725 726 #undef __FUNCT__ 727 #define __FUNCT__ "SNESTSFormJacobian_RosW" 728 static PetscErrorCode SNESTSFormJacobian_RosW(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *str,TS ts) 729 { 730 TS_RosW *ros = (TS_RosW*)ts->data; 731 PetscErrorCode ierr; 732 733 PetscFunctionBegin; 734 /* ros->Ydot and ros->Ystage have already been computed in SNESTSFormFunction_RosW (SNES guarantees this) */ 735 ierr = TSComputeIJacobian(ts,ros->stage_time,ros->Ystage,ros->Ydot,ros->shift,A,B,str,PETSC_TRUE);CHKERRQ(ierr); 736 PetscFunctionReturn(0); 737 } 738 739 #undef __FUNCT__ 740 #define __FUNCT__ "TSSetUp_RosW" 741 static PetscErrorCode TSSetUp_RosW(TS ts) 742 { 743 TS_RosW *ros = (TS_RosW*)ts->data; 744 RosWTableau tab = ros->tableau; 745 PetscInt s = tab->s; 746 PetscErrorCode ierr; 747 748 PetscFunctionBegin; 749 if (!ros->tableau) { 750 ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr); 751 } 752 ierr = VecDuplicateVecs(ts->vec_sol,s,&ros->Y);CHKERRQ(ierr); 753 ierr = VecDuplicate(ts->vec_sol,&ros->Ydot);CHKERRQ(ierr); 754 ierr = VecDuplicate(ts->vec_sol,&ros->Ystage);CHKERRQ(ierr); 755 ierr = VecDuplicate(ts->vec_sol,&ros->Zdot);CHKERRQ(ierr); 756 ierr = VecDuplicate(ts->vec_sol,&ros->Zstage);CHKERRQ(ierr); 757 ierr = PetscMalloc(s*sizeof(ros->work[0]),&ros->work);CHKERRQ(ierr); 758 PetscFunctionReturn(0); 759 } 760 /*------------------------------------------------------------*/ 761 762 #undef __FUNCT__ 763 #define __FUNCT__ "TSSetFromOptions_RosW" 764 static PetscErrorCode TSSetFromOptions_RosW(TS ts) 765 { 766 TS_RosW *ros = (TS_RosW*)ts->data; 767 PetscErrorCode ierr; 768 char rostype[256]; 769 770 PetscFunctionBegin; 771 ierr = PetscOptionsHead("RosW ODE solver options");CHKERRQ(ierr); 772 { 773 RosWTableauLink link; 774 PetscInt count,choice; 775 PetscBool flg; 776 const char **namelist; 777 SNES snes; 778 779 ierr = PetscStrncpy(rostype,TSRosWDefault,sizeof rostype);CHKERRQ(ierr); 780 for (link=RosWTableauList,count=0; link; link=link->next,count++) ; 781 ierr = PetscMalloc(count*sizeof(char*),&namelist);CHKERRQ(ierr); 782 for (link=RosWTableauList,count=0; link; link=link->next,count++) namelist[count] = link->tab.name; 783 ierr = PetscOptionsEList("-ts_rosw_type","Family of Rosenbrock-W method","TSRosWSetType",(const char*const*)namelist,count,rostype,&choice,&flg);CHKERRQ(ierr); 784 ierr = TSRosWSetType(ts,flg ? namelist[choice] : rostype);CHKERRQ(ierr); 785 ierr = PetscFree(namelist);CHKERRQ(ierr); 786 787 ierr = PetscOptionsBool("-ts_rosw_recompute_jacobian","Recompute the Jacobian at each stage","TSRosWSetRecomputeJacobian",ros->recompute_jacobian,&ros->recompute_jacobian,PETSC_NULL);CHKERRQ(ierr); 788 789 /* Rosenbrock methods are linearly implicit, so set that unless the user has specifically asked for something else */ 790 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 791 if (!((PetscObject)snes)->type_name) { 792 ierr = SNESSetType(snes,SNESKSPONLY);CHKERRQ(ierr); 793 } 794 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 795 } 796 ierr = PetscOptionsTail();CHKERRQ(ierr); 797 PetscFunctionReturn(0); 798 } 799 800 #undef __FUNCT__ 801 #define __FUNCT__ "PetscFormatRealArray" 802 static PetscErrorCode PetscFormatRealArray(char buf[],size_t len,const char *fmt,PetscInt n,const PetscReal x[]) 803 { 804 PetscErrorCode ierr; 805 PetscInt i; 806 size_t left,count; 807 char *p; 808 809 PetscFunctionBegin; 810 for (i=0,p=buf,left=len; i<n; i++) { 811 ierr = PetscSNPrintfCount(p,left,fmt,&count,x[i]);CHKERRQ(ierr); 812 if (count >= left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Insufficient space in buffer"); 813 left -= count; 814 p += count; 815 *p++ = ' '; 816 } 817 p[i ? 0 : -1] = 0; 818 PetscFunctionReturn(0); 819 } 820 821 #undef __FUNCT__ 822 #define __FUNCT__ "TSView_RosW" 823 static PetscErrorCode TSView_RosW(TS ts,PetscViewer viewer) 824 { 825 TS_RosW *ros = (TS_RosW*)ts->data; 826 RosWTableau tab = ros->tableau; 827 PetscBool iascii; 828 PetscErrorCode ierr; 829 830 PetscFunctionBegin; 831 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 832 if (iascii) { 833 const TSRosWType rostype; 834 PetscInt i; 835 PetscReal abscissa[512]; 836 char buf[512]; 837 ierr = TSRosWGetType(ts,&rostype);CHKERRQ(ierr); 838 ierr = PetscViewerASCIIPrintf(viewer," Rosenbrock-W %s\n",rostype);CHKERRQ(ierr); 839 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,tab->ASum);CHKERRQ(ierr); 840 ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A = %s\n",buf);CHKERRQ(ierr); 841 for (i=0; i<tab->s; i++) abscissa[i] = tab->ASum[i] + tab->Gamma[i]; 842 ierr = PetscFormatRealArray(buf,sizeof buf,"% 8.6f",tab->s,abscissa);CHKERRQ(ierr); 843 ierr = PetscViewerASCIIPrintf(viewer," Abscissa of A+Gamma = %s\n",buf);CHKERRQ(ierr); 844 } 845 ierr = SNESView(ts->snes,viewer);CHKERRQ(ierr); 846 PetscFunctionReturn(0); 847 } 848 849 #undef __FUNCT__ 850 #define __FUNCT__ "TSRosWSetType" 851 /*@C 852 TSRosWSetType - Set the type of Rosenbrock-W scheme 853 854 Logically collective 855 856 Input Parameter: 857 + ts - timestepping context 858 - rostype - type of Rosenbrock-W scheme 859 860 Level: intermediate 861 862 .seealso: TSRosWGetType() 863 @*/ 864 PetscErrorCode TSRosWSetType(TS ts,const TSRosWType rostype) 865 { 866 PetscErrorCode ierr; 867 868 PetscFunctionBegin; 869 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 870 ierr = PetscTryMethod(ts,"TSRosWSetType_C",(TS,const TSRosWType),(ts,rostype));CHKERRQ(ierr); 871 PetscFunctionReturn(0); 872 } 873 874 #undef __FUNCT__ 875 #define __FUNCT__ "TSRosWGetType" 876 /*@C 877 TSRosWGetType - Get the type of Rosenbrock-W scheme 878 879 Logically collective 880 881 Input Parameter: 882 . ts - timestepping context 883 884 Output Parameter: 885 . rostype - type of Rosenbrock-W scheme 886 887 Level: intermediate 888 889 .seealso: TSRosWGetType() 890 @*/ 891 PetscErrorCode TSRosWGetType(TS ts,const TSRosWType *rostype) 892 { 893 PetscErrorCode ierr; 894 895 PetscFunctionBegin; 896 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 897 ierr = PetscUseMethod(ts,"TSRosWGetType_C",(TS,const TSRosWType*),(ts,rostype));CHKERRQ(ierr); 898 PetscFunctionReturn(0); 899 } 900 901 #undef __FUNCT__ 902 #define __FUNCT__ "TSRosWSetRecomputeJacobian" 903 /*@C 904 TSRosWSetRecomputeJacobian - Set whether to recompute the Jacobian at each stage. The default is to update the Jacobian once per step. 905 906 Logically collective 907 908 Input Parameter: 909 + ts - timestepping context 910 - flg - PETSC_TRUE to recompute the Jacobian at each stage 911 912 Level: intermediate 913 914 .seealso: TSRosWGetType() 915 @*/ 916 PetscErrorCode TSRosWSetRecomputeJacobian(TS ts,PetscBool flg) 917 { 918 PetscErrorCode ierr; 919 920 PetscFunctionBegin; 921 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 922 ierr = PetscTryMethod(ts,"TSRosWSetRecomputeJacobian_C",(TS,PetscBool),(ts,flg));CHKERRQ(ierr); 923 PetscFunctionReturn(0); 924 } 925 926 EXTERN_C_BEGIN 927 #undef __FUNCT__ 928 #define __FUNCT__ "TSRosWGetType_RosW" 929 PetscErrorCode TSRosWGetType_RosW(TS ts,const TSRosWType *rostype) 930 { 931 TS_RosW *ros = (TS_RosW*)ts->data; 932 PetscErrorCode ierr; 933 934 PetscFunctionBegin; 935 if (!ros->tableau) {ierr = TSRosWSetType(ts,TSRosWDefault);CHKERRQ(ierr);} 936 *rostype = ros->tableau->name; 937 PetscFunctionReturn(0); 938 } 939 #undef __FUNCT__ 940 #define __FUNCT__ "TSRosWSetType_RosW" 941 PetscErrorCode TSRosWSetType_RosW(TS ts,const TSRosWType rostype) 942 { 943 TS_RosW *ros = (TS_RosW*)ts->data; 944 PetscErrorCode ierr; 945 PetscBool match; 946 RosWTableauLink link; 947 948 PetscFunctionBegin; 949 if (ros->tableau) { 950 ierr = PetscStrcmp(ros->tableau->name,rostype,&match);CHKERRQ(ierr); 951 if (match) PetscFunctionReturn(0); 952 } 953 for (link = RosWTableauList; link; link=link->next) { 954 ierr = PetscStrcmp(link->tab.name,rostype,&match);CHKERRQ(ierr); 955 if (match) { 956 ierr = TSReset_RosW(ts);CHKERRQ(ierr); 957 ros->tableau = &link->tab; 958 PetscFunctionReturn(0); 959 } 960 } 961 SETERRQ1(((PetscObject)ts)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rostype); 962 PetscFunctionReturn(0); 963 } 964 965 #undef __FUNCT__ 966 #define __FUNCT__ "TSRosWSetRecomputeJacobian_RosW" 967 PetscErrorCode TSRosWSetRecomputeJacobian_RosW(TS ts,PetscBool flg) 968 { 969 TS_RosW *ros = (TS_RosW*)ts->data; 970 971 PetscFunctionBegin; 972 ros->recompute_jacobian = flg; 973 PetscFunctionReturn(0); 974 } 975 EXTERN_C_END 976 977 /* ------------------------------------------------------------ */ 978 /*MC 979 TSRosW - ODE solver using Rosenbrock-W schemes 980 981 These methods are intended for problems with well-separated time scales, especially when a slow scale is strongly 982 nonlinear such that it is expensive to solve with a fully implicit method. The user should provide the stiff part 983 of the equation using TSSetIFunction() and the non-stiff part with TSSetRHSFunction(). 984 985 Notes: 986 This method currently only works with autonomous ODE and DAE. 987 988 Developer notes: 989 Rosenbrock-W methods are typically specified for autonomous ODE 990 991 $ xdot = f(x) 992 993 by the stage equations 994 995 $ k_i = h f(x_0 + sum_j alpha_ij k_j) + h J sum_j gamma_ij k_j 996 997 and step completion formula 998 999 $ x_1 = x_0 + sum_j b_j k_j 1000 1001 with step size h and coefficients alpha_ij, gamma_ij, and b_i. Implementing the method in this form would require f(x) 1002 and the Jacobian J to be available, in addition to the shifted matrix I - h gamma_ii J. Following Hairer and Wanner, 1003 we define new variables for the stage equations 1004 1005 $ y_i = gamma_ij k_j 1006 1007 The k_j can be recovered because Gamma is invertible. Let C be the lower triangular part of Gamma^{-1} and define 1008 1009 $ A = Alpha Gamma^{-1}, bt^T = b^T Gamma^{-i} 1010 1011 to rewrite the method as 1012 1013 $ [M/(h gamma_ii) - J] y_i = f(x_0 + sum_j a_ij y_j) + M sum_j (c_ij/h) y_j 1014 $ x_1 = x_0 + sum_j bt_j y_j 1015 1016 where we have introduced the mass matrix M. Continue by defining 1017 1018 $ ydot_i = 1/(h gamma_ii) y_i - sum_j (c_ij/h) y_j 1019 1020 or, more compactly in tensor notation 1021 1022 $ Ydot = 1/h (Gamma^{-1} \otimes I) Y . 1023 1024 Note that Gamma^{-1} is lower triangular. With this definition of Ydot in terms of known quantities and the current 1025 stage y_i, the stage equations reduce to performing one Newton step (typically with a lagged Jacobian) on the 1026 equation 1027 1028 $ g(x_0 + sum_j a_ij y_j + y_i, ydot_i) = 0 1029 1030 with initial guess y_i = 0. 1031 1032 Level: beginner 1033 1034 .seealso: TSCreate(), TS, TSSetType(), TSRosWRegister() 1035 1036 M*/ 1037 EXTERN_C_BEGIN 1038 #undef __FUNCT__ 1039 #define __FUNCT__ "TSCreate_RosW" 1040 PetscErrorCode TSCreate_RosW(TS ts) 1041 { 1042 TS_RosW *ros; 1043 PetscErrorCode ierr; 1044 1045 PetscFunctionBegin; 1046 #if !defined(PETSC_USE_DYNAMIC_LIBRARIES) 1047 ierr = TSRosWInitializePackage(PETSC_NULL);CHKERRQ(ierr); 1048 #endif 1049 1050 ts->ops->reset = TSReset_RosW; 1051 ts->ops->destroy = TSDestroy_RosW; 1052 ts->ops->view = TSView_RosW; 1053 ts->ops->setup = TSSetUp_RosW; 1054 ts->ops->step = TSStep_RosW; 1055 ts->ops->interpolate = TSInterpolate_RosW; 1056 ts->ops->evaluatestep = TSEvaluateStep_RosW; 1057 ts->ops->setfromoptions = TSSetFromOptions_RosW; 1058 ts->ops->snesfunction = SNESTSFormFunction_RosW; 1059 ts->ops->snesjacobian = SNESTSFormJacobian_RosW; 1060 1061 ierr = PetscNewLog(ts,TS_RosW,&ros);CHKERRQ(ierr); 1062 ts->data = (void*)ros; 1063 1064 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWGetType_C","TSRosWGetType_RosW",TSRosWGetType_RosW);CHKERRQ(ierr); 1065 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetType_C","TSRosWSetType_RosW",TSRosWSetType_RosW);CHKERRQ(ierr); 1066 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRosWSetRecomputeJacobian_C","TSRosWSetRecomputeJacobian_RosW",TSRosWSetRecomputeJacobian_RosW);CHKERRQ(ierr); 1067 PetscFunctionReturn(0); 1068 } 1069 EXTERN_C_END 1070