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