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