Actual source code: ex9.c

  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: {
 19:   VecWAXPY(F, 1.0, U, V);
 20:   return 0;
 21: }

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

 26: */
 27: PetscErrorCode F(PetscReal t, Vec U, Vec V, Vec F)
 28: {
 30:   VecWAXPY(F, -1.0, V, U);
 31:   return 0;
 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;

 55:   PetscInitialize(&argc, &argv, (char *)0, help);
 56:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 57:   TSCreate(PETSC_COMM_WORLD, &ts);
 58:   TSSetProblemType(ts, TS_NONLINEAR);
 59:   TSSetType(ts, TSROSW);
 60:   TSSetFromOptions(ts);
 61:   VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &tsrhs);
 62:   VecCreateMPI(PETSC_COMM_WORLD, 2, PETSC_DETERMINE, &UV);
 63:   TSSetRHSFunction(ts, tsrhs, TSFunctionRHS, &ctx);
 64:   TSSetIFunction(ts, NULL, TSFunctionI, &ctx);
 65:   ctx.f = f;
 66:   ctx.F = F;

 68:   VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.U);
 69:   VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.V);
 70:   VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.UF);
 71:   VecCreateMPI(PETSC_COMM_WORLD, 1, PETSC_DETERMINE, &ctx.VF);
 72:   I = 2 * rank;
 73:   ISCreateGeneral(PETSC_COMM_WORLD, 1, &I, PETSC_COPY_VALUES, &is);
 74:   VecScatterCreate(ctx.U, NULL, UV, is, &ctx.scatterU);
 75:   ISDestroy(&is);
 76:   I = 2 * rank + 1;
 77:   ISCreateGeneral(PETSC_COMM_WORLD, 1, &I, PETSC_COPY_VALUES, &is);
 78:   VecScatterCreate(ctx.V, NULL, UV, is, &ctx.scatterV);
 79:   ISDestroy(&is);

 81:   VecSet(UV, 1.0);
 82:   TSSolve(ts, UV);
 83:   VecDestroy(&tsrhs);
 84:   VecDestroy(&UV);
 85:   VecDestroy(&ctx.U);
 86:   VecDestroy(&ctx.V);
 87:   VecDestroy(&ctx.UF);
 88:   VecDestroy(&ctx.VF);
 89:   VecScatterDestroy(&ctx.scatterU);
 90:   VecScatterDestroy(&ctx.scatterV);
 91:   TSDestroy(&ts);
 92:   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;

105:   VecSet(F, 0.0);
106:   VecScatterBegin(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE);
107:   VecScatterEnd(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE);
108:   VecScatterBegin(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE);
109:   VecScatterEnd(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE);
110:   (*ctx->f)(t, ctx->U, ctx->V, ctx->UF);
111:   VecScatterBegin(ctx->scatterU, ctx->UF, F, INSERT_VALUES, SCATTER_FORWARD);
112:   VecScatterEnd(ctx->scatterU, ctx->UF, F, INSERT_VALUES, SCATTER_FORWARD);
113:   return 0;
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;

125:   VecCopy(UVdot, F);
126:   VecScatterBegin(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE);
127:   VecScatterEnd(ctx->scatterU, UV, ctx->U, INSERT_VALUES, SCATTER_REVERSE);
128:   VecScatterBegin(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE);
129:   VecScatterEnd(ctx->scatterV, UV, ctx->V, INSERT_VALUES, SCATTER_REVERSE);
130:   (*ctx->F)(t, ctx->U, ctx->V, ctx->VF);
131:   VecScatterBegin(ctx->scatterV, ctx->VF, F, INSERT_VALUES, SCATTER_FORWARD);
132:   VecScatterEnd(ctx->scatterV, ctx->VF, F, INSERT_VALUES, SCATTER_FORWARD);
133:   return 0;
134: }