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 Same as ex6.c and ex7.c except calls the ARKIMEX integrator on the entire DAE 10 */ 11 12 /* 13 f(U,V) = U + V 14 15 */ 16 PetscErrorCode f(PetscReal t, Vec U, Vec V, Vec F) { 17 PetscFunctionBeginUser; 18 PetscCall(VecWAXPY(F, 1.0, U, V)); 19 PetscFunctionReturn(0); 20 } 21 22 /* 23 F(U,V) = U - V 24 25 */ 26 PetscErrorCode F(PetscReal t, Vec U, Vec V, Vec F) { 27 PetscFunctionBeginUser; 28 PetscCall(VecWAXPY(F, -1.0, V, U)); 29 PetscFunctionReturn(0); 30 } 31 32 typedef struct { 33 Vec U, V; 34 Vec UF, VF; 35 VecScatter scatterU, scatterV; 36 PetscErrorCode (*f)(PetscReal, Vec, Vec, Vec); 37 PetscErrorCode (*F)(PetscReal, Vec, Vec, Vec); 38 } AppCtx; 39 40 extern PetscErrorCode TSFunctionRHS(TS, PetscReal, Vec, Vec, void *); 41 extern PetscErrorCode TSFunctionI(TS, PetscReal, Vec, Vec, Vec, void *); 42 43 int main(int argc, char **argv) { 44 AppCtx ctx; 45 TS ts; 46 Vec tsrhs, UV; 47 IS is; 48 PetscInt I; 49 PetscMPIInt rank; 50 51 PetscFunctionBeginUser; 52 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 53 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 54 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 55 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 56 PetscCall(TSSetType(ts, TSROSW)); 57 PetscCall(TSSetFromOptions(ts)); 58 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &tsrhs)); 59 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &UV)); 60 PetscCall(TSSetRHSFunction(ts, tsrhs, TSFunctionRHS, &ctx)); 61 PetscCall(TSSetIFunction(ts, NULL, TSFunctionI, &ctx)); 62 ctx.f = f; 63 ctx.F = F; 64 65 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.U)); 66 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.V)); 67 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.UF)); 68 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.VF)); 69 I = 2 * rank; 70 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 1, &I, PETSC_COPY_VALUES, &is)); 71 PetscCall(VecScatterCreate(ctx.U, NULL, UV, is, &ctx.scatterU)); 72 PetscCall(ISDestroy(&is)); 73 I = 2 * rank + 1; 74 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 1, &I, PETSC_COPY_VALUES, &is)); 75 PetscCall(VecScatterCreate(ctx.V, NULL, UV, is, &ctx.scatterV)); 76 PetscCall(ISDestroy(&is)); 77 78 PetscCall(VecSet(UV, 1.0)); 79 PetscCall(TSSolve(ts, UV)); 80 PetscCall(VecDestroy(&tsrhs)); 81 PetscCall(VecDestroy(&UV)); 82 PetscCall(VecDestroy(&ctx.U)); 83 PetscCall(VecDestroy(&ctx.V)); 84 PetscCall(VecDestroy(&ctx.UF)); 85 PetscCall(VecDestroy(&ctx.VF)); 86 PetscCall(VecScatterDestroy(&ctx.scatterU)); 87 PetscCall(VecScatterDestroy(&ctx.scatterV)); 88 PetscCall(TSDestroy(&ts)); 89 PetscCall(PetscFinalize()); 90 return 0; 91 } 92 93 /* 94 Defines the RHS function that is passed to the time-integrator. 95 96 */ 97 PetscErrorCode TSFunctionRHS(TS ts, PetscReal t, Vec UV, Vec F, void *actx) { 98 AppCtx *ctx = (AppCtx *)actx; 99 100 PetscFunctionBeginUser; 101 PetscCall(VecSet(F, 0.0)); 102 PetscCall(VecScatterBegin(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE)); 103 PetscCall(VecScatterEnd(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE)); 104 PetscCall(VecScatterBegin(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE)); 105 PetscCall(VecScatterEnd(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE)); 106 PetscCall((*ctx->f)(t, ctx->U, ctx->V, ctx->UF)); 107 PetscCall(VecScatterBegin(ctx->scatterU, ctx->UF, F, INSERT_VALUES, SCATTER_FORWARD)); 108 PetscCall(VecScatterEnd(ctx->scatterU, ctx->UF, F, INSERT_VALUES, SCATTER_FORWARD)); 109 PetscFunctionReturn(0); 110 } 111 112 /* 113 Defines the nonlinear function that is passed to the time-integrator 114 115 */ 116 PetscErrorCode TSFunctionI(TS ts, PetscReal t, Vec UV, Vec UVdot, Vec F, void *actx) { 117 AppCtx *ctx = (AppCtx *)actx; 118 119 PetscFunctionBeginUser; 120 PetscCall(VecCopy(UVdot, F)); 121 PetscCall(VecScatterBegin(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE)); 122 PetscCall(VecScatterEnd(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE)); 123 PetscCall(VecScatterBegin(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE)); 124 PetscCall(VecScatterEnd(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE)); 125 PetscCall((*ctx->F)(t, ctx->U, ctx->V, ctx->VF)); 126 PetscCall(VecScatterBegin(ctx->scatterV, ctx->VF, F, INSERT_VALUES, SCATTER_FORWARD)); 127 PetscCall(VecScatterEnd(ctx->scatterV, ctx->VF, F, INSERT_VALUES, SCATTER_FORWARD)); 128 PetscFunctionReturn(0); 129 } 130