1 static char help[] = "Solves DAE with integrator only on non-algebraic terms \n"; 2 3 #include <petscts.h> 4 5 /* 6 \dot{U} = f(U,V) 7 F(U,V) = 0 8 9 Same as ex6.c and ex7.c except calls the TSROSW integrator on the entire DAE 10 */ 11 12 /* 13 f(U,V) = U + V 14 15 */ 16 PetscErrorCode f(PetscReal t,Vec UV,Vec F) 17 { 18 const PetscScalar *u,*v; 19 PetscScalar *f; 20 PetscInt n,i; 21 22 PetscFunctionBeginUser; 23 CHKERRQ(VecGetLocalSize(UV,&n)); 24 n = n/2; 25 CHKERRQ(VecGetArrayRead(UV,&u)); 26 v = u + n; 27 CHKERRQ(VecGetArrayWrite(F,&f)); 28 for (i=0; i<n; i++) f[i] = u[i] + v[i]; 29 CHKERRQ(VecRestoreArrayRead(UV,&u)); 30 CHKERRQ(VecRestoreArrayWrite(F,&f)); 31 PetscFunctionReturn(0); 32 } 33 34 /* 35 F(U,V) = U - V 36 37 */ 38 PetscErrorCode F(PetscReal t,Vec UV,Vec F) 39 { 40 const PetscScalar *u,*v; 41 PetscScalar *f; 42 PetscInt n,i; 43 44 PetscFunctionBeginUser; 45 CHKERRQ(VecGetLocalSize(UV,&n)); 46 n = n/2; 47 CHKERRQ(VecGetArrayRead(UV,&u)); 48 v = u + n; 49 CHKERRQ(VecGetArrayWrite(F,&f)); 50 f = f + n; 51 for (i=0; i<n; i++) f[i] = u[i] - v[i]; 52 f = f - n; 53 CHKERRQ(VecRestoreArrayRead(UV,&u)); 54 CHKERRQ(VecRestoreArrayWrite(F,&f)); 55 PetscFunctionReturn(0); 56 } 57 58 typedef struct { 59 PetscErrorCode (*f)(PetscReal,Vec,Vec); 60 PetscErrorCode (*F)(PetscReal,Vec,Vec); 61 } AppCtx; 62 63 extern PetscErrorCode TSFunctionRHS(TS,PetscReal,Vec,Vec,void*); 64 extern PetscErrorCode TSFunctionI(TS,PetscReal,Vec,Vec,Vec,void*); 65 66 int main(int argc,char **argv) 67 { 68 AppCtx ctx; 69 TS ts; 70 Vec tsrhs,UV; 71 72 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 73 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 74 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 75 CHKERRQ(TSSetType(ts,TSROSW)); 76 CHKERRQ(TSSetFromOptions(ts)); 77 CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&tsrhs)); 78 CHKERRQ(VecDuplicate(tsrhs,&UV)); 79 CHKERRQ(TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx)); 80 CHKERRQ(TSSetIFunction(ts,NULL,TSFunctionI,&ctx)); 81 CHKERRQ(TSSetMaxTime(ts,1.0)); 82 ctx.f = f; 83 ctx.F = F; 84 85 CHKERRQ(VecSet(UV,1.0)); 86 CHKERRQ(TSSolve(ts,UV)); 87 CHKERRQ(VecDestroy(&tsrhs)); 88 CHKERRQ(VecDestroy(&UV)); 89 CHKERRQ(TSDestroy(&ts)); 90 CHKERRQ(PetscFinalize()); 91 return 0; 92 } 93 94 /* 95 Defines the RHS function that is passed to the time-integrator. 96 */ 97 PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx) 98 { 99 AppCtx *ctx = (AppCtx*)actx; 100 101 PetscFunctionBeginUser; 102 CHKERRQ(VecSet(F,0.0)); 103 CHKERRQ((*ctx->f)(t,UV,F)); 104 PetscFunctionReturn(0); 105 } 106 107 /* 108 Defines the nonlinear function that is passed to the time-integrator 109 */ 110 PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx) 111 { 112 AppCtx *ctx = (AppCtx*)actx; 113 114 PetscFunctionBeginUser; 115 CHKERRQ(VecCopy(UVdot,F)); 116 CHKERRQ((*ctx->F)(t,UV,F)); 117 PetscFunctionReturn(0); 118 } 119 120 /*TEST 121 122 test: 123 args: -ts_view 124 125 test: 126 suffix: 2 127 args: -snes_lag_jacobian 2 -ts_view 128 129 TEST*/ 130