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