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