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 PetscCall(VecWAXPY(F, 1.0, U, V)); 17 PetscFunctionReturn(PETSC_SUCCESS); 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 PetscCall(VecWAXPY(F, -1.0, V, U)); 27 PetscFunctionReturn(PETSC_SUCCESS); 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 AppCtx ctx; 44 TS ts; 45 Vec tsrhs, U; 46 47 PetscFunctionBeginUser; 48 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 49 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 50 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 51 PetscCall(TSSetType(ts, TSEULER)); 52 PetscCall(TSSetFromOptions(ts)); 53 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 1, &tsrhs)); 54 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 1, &U)); 55 PetscCall(TSSetRHSFunction(ts, tsrhs, TSFunction, &ctx)); 56 PetscCall(TSSetMaxTime(ts, 1.0)); 57 ctx.f = f; 58 59 PetscCall(SNESCreate(PETSC_COMM_WORLD, &ctx.snes)); 60 PetscCall(SNESSetFromOptions(ctx.snes)); 61 PetscCall(SNESSetFunction(ctx.snes, NULL, SNESFunction, &ctx)); 62 PetscCall(SNESSetJacobian(ctx.snes, NULL, NULL, SNESComputeJacobianDefault, &ctx)); 63 ctx.F = F; 64 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 1, &ctx.V)); 65 66 PetscCall(VecSet(U, 1.0)); 67 PetscCall(TSSolve(ts, U)); 68 69 PetscCall(VecDestroy(&ctx.V)); 70 PetscCall(VecDestroy(&tsrhs)); 71 PetscCall(VecDestroy(&U)); 72 PetscCall(SNESDestroy(&ctx.snes)); 73 PetscCall(TSDestroy(&ts)); 74 PetscCall(PetscFinalize()); 75 return 0; 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 PetscCall(SNESSolve(ctx->snes, NULL, ctx->V)); 91 PetscCall((*ctx->f)(t, U, ctx->V, F)); 92 PetscFunctionReturn(PETSC_SUCCESS); 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 PetscCall((*ctx->F)(ctx->t, ctx->U, V, F)); 104 PetscFunctionReturn(PETSC_SUCCESS); 105 } 106 107 /*TEST 108 109 test: 110 args: -ts_monitor -ts_view 111 112 TEST*/ 113