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 10 /* 11 f(U,V) = U + V 12 */ 13 PetscErrorCode f(PetscReal t,Vec U,Vec V,Vec F) 14 { 15 PetscFunctionBeginUser; 16 CHKERRQ(VecWAXPY(F,1.0,U,V)); 17 PetscFunctionReturn(0); 18 } 19 20 /* 21 F(U,V) = U - V 22 */ 23 PetscErrorCode F(PetscReal t,Vec U,Vec V,Vec F) 24 { 25 PetscFunctionBeginUser; 26 CHKERRQ(VecWAXPY(F,-1.0,V,U)); 27 PetscFunctionReturn(0); 28 } 29 30 typedef struct { 31 PetscReal t; 32 SNES snes; 33 Vec U,V; 34 PetscErrorCode (*f)(PetscReal,Vec,Vec,Vec); 35 PetscErrorCode (*F)(PetscReal,Vec,Vec,Vec); 36 } AppCtx; 37 38 extern PetscErrorCode TSFunction(TS,PetscReal,Vec,Vec,void*); 39 extern PetscErrorCode SNESFunction(SNES,Vec,Vec,void*); 40 41 int main(int argc,char **argv) 42 { 43 PetscErrorCode ierr; 44 AppCtx ctx; 45 TS ts; 46 Vec tsrhs,U; 47 48 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 49 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 50 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 51 CHKERRQ(TSSetType(ts,TSEULER)); 52 CHKERRQ(TSSetFromOptions(ts)); 53 CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&tsrhs)); 54 CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&U)); 55 CHKERRQ(TSSetRHSFunction(ts,tsrhs,TSFunction,&ctx)); 56 CHKERRQ(TSSetMaxTime(ts,1.0)); 57 ctx.f = f; 58 59 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&ctx.snes)); 60 CHKERRQ(SNESSetFromOptions(ctx.snes)); 61 CHKERRQ(SNESSetFunction(ctx.snes,NULL,SNESFunction,&ctx)); 62 CHKERRQ(SNESSetJacobian(ctx.snes,NULL,NULL,SNESComputeJacobianDefault,&ctx)); 63 ctx.F = F; 64 CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,1,&ctx.V)); 65 66 CHKERRQ(VecSet(U,1.0)); 67 CHKERRQ(TSSolve(ts,U)); 68 69 CHKERRQ(VecDestroy(&ctx.V)); 70 CHKERRQ(VecDestroy(&tsrhs)); 71 CHKERRQ(VecDestroy(&U)); 72 CHKERRQ(SNESDestroy(&ctx.snes)); 73 CHKERRQ(TSDestroy(&ts)); 74 ierr = PetscFinalize(); 75 return ierr; 76 } 77 78 /* 79 Defines the RHS function that is passed to the time-integrator. 80 81 Solves F(U,V) for V and then computes f(U,V) 82 */ 83 PetscErrorCode TSFunction(TS ts,PetscReal t,Vec U,Vec F,void *actx) 84 { 85 AppCtx *ctx = (AppCtx*)actx; 86 87 PetscFunctionBeginUser; 88 ctx->t = t; 89 ctx->U = U; 90 CHKERRQ(SNESSolve(ctx->snes,NULL,ctx->V)); 91 CHKERRQ((*ctx->f)(t,U,ctx->V,F)); 92 PetscFunctionReturn(0); 93 } 94 95 /* 96 Defines the nonlinear function that is passed to the nonlinear solver 97 */ 98 PetscErrorCode SNESFunction(SNES snes,Vec V,Vec F,void *actx) 99 { 100 AppCtx *ctx = (AppCtx*)actx; 101 102 PetscFunctionBeginUser; 103 CHKERRQ((*ctx->F)(ctx->t,ctx->U,V,F)); 104 PetscFunctionReturn(0); 105 } 106 107 /*TEST 108 109 test: 110 args: -ts_monitor -ts_view 111 112 TEST*/ 113