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