1: static char help[] = "Solves DAE with integrator only on non-algebraic terms \n";

  3: #include <petscts.h>

  5: /*
  6:         \dot{U} = f(U,V)
  7:         F(U,V)  = 0

  9:     Same as ex6.c and ex7.c except calls the ARKIMEX integrator on the entire DAE
 10: */

 12: /*
 13:    f(U,V) = U + V

 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: }

 23: /*
 24:    F(U,V) = U - V

 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: }

 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;

 42: extern PetscErrorCode TSFunctionRHS(TS, PetscReal, Vec, Vec, void *);
 43: extern PetscErrorCode TSFunctionI(TS, PetscReal, Vec, Vec, Vec, void *);

 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;

 54:   PetscFunctionBeginUser;
 55:   PetscCall(PetscInitialize(&argc, &argv, NULL, 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;

 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));

 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: }

 96: /*
 97:    Defines the RHS function that is passed to the time-integrator.

 99: */
100: PetscErrorCode TSFunctionRHS(TS ts, PetscReal t, Vec UV, Vec F, void *actx)
101: {
102:   AppCtx *ctx = (AppCtx *)actx;

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: }

116: /*
117:    Defines the nonlinear function that is passed to the time-integrator

119: */
120: PetscErrorCode TSFunctionI(TS ts, PetscReal t, Vec UV, Vec UVdot, Vec F, void *actx)
121: {
122:   AppCtx *ctx = (AppCtx *)actx;

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: }