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