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