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 ARKIMEX 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 = VecGetArray(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 = VecRestoreArray(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 = VecGetArray(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 = VecRestoreArray(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 = VecCreateMPI(PETSC_COMM_WORLD,2,PETSC_DETERMINE,&UV);CHKERRQ(ierr); 83 ierr = TSSetRHSFunction(ts,tsrhs,TSFunctionRHS,&ctx);CHKERRQ(ierr); 84 ierr = TSSetIFunction(ts,NULL,TSFunctionI,&ctx);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 101 102 */ 103 PetscErrorCode TSFunctionRHS(TS ts,PetscReal t,Vec UV,Vec F,void *actx) 104 { 105 AppCtx *ctx = (AppCtx*)actx; 106 PetscErrorCode ierr; 107 108 PetscFunctionBeginUser; 109 ierr = VecSet(F,0.0);CHKERRQ(ierr); 110 ierr = (*ctx->f)(t,UV,F);CHKERRQ(ierr); 111 PetscFunctionReturn(0); 112 } 113 114 /* 115 Defines the nonlinear function that is passed to the time-integrator 116 117 */ 118 PetscErrorCode TSFunctionI(TS ts,PetscReal t,Vec UV,Vec UVdot,Vec F,void *actx) 119 { 120 AppCtx *ctx = (AppCtx*)actx; 121 PetscErrorCode ierr; 122 123 PetscFunctionBeginUser; 124 ierr = VecCopy(UVdot,F);CHKERRQ(ierr); 125 ierr = (*ctx->F)(t,UV,F);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 129 130