Actual source code: ex1.c

  1: static char help[] = "Basic problem for multi-rate method.\n";

  3: /*F

  5: \begin{eqnarray}
  6:                  ys' = \frac{ys}{a}\\
  7:                  yf' = ys*cos(bt)\\
  8: \end{eqnarray}

 10: F*/

 12: #include <petscts.h>

 14: typedef struct {
 15:   PetscReal a, b, Tf, dt;
 16: } AppCtx;

 18: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
 19: {
 20:   const PetscScalar *u;
 21:   PetscScalar       *f;

 23:   VecGetArrayRead(U, &u);
 24:   VecGetArray(F, &f);
 25:   f[0] = u[0] / ctx->a;
 26:   f[1] = u[0] * PetscCosScalar(t * ctx->b);
 27:   VecRestoreArrayRead(U, &u);
 28:   VecRestoreArray(F, &f);
 29:   return 0;
 30: }

 32: static PetscErrorCode RHSFunctionslow(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
 33: {
 34:   const PetscScalar *u;
 35:   PetscScalar       *f;

 37:   VecGetArrayRead(U, &u);
 38:   VecGetArray(F, &f);
 39:   f[0] = u[0] / ctx->a;
 40:   VecRestoreArrayRead(U, &u);
 41:   VecRestoreArray(F, &f);
 42:   return 0;
 43: }

 45: static PetscErrorCode RHSFunctionfast(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
 46: {
 47:   const PetscScalar *u;
 48:   PetscScalar       *f;

 50:   VecGetArrayRead(U, &u);
 51:   VecGetArray(F, &f);
 52:   f[0] = u[0] * PetscCosScalar(t * ctx->b);
 53:   VecRestoreArrayRead(U, &u);
 54:   VecRestoreArray(F, &f);
 55:   return 0;
 56: }

 58: /*
 59:   Define the analytic solution for check method easily
 60: */
 61: static PetscErrorCode sol_true(PetscReal t, Vec U, AppCtx *ctx)
 62: {
 63:   PetscScalar *u;

 65:   VecGetArray(U, &u);
 66:   u[0] = PetscExpScalar(t / ctx->a);
 67:   u[1] = (ctx->a * PetscCosScalar(ctx->b * t) + ctx->a * ctx->a * ctx->b * PetscSinScalar(ctx->b * t)) * PetscExpScalar(t / ctx->a) / (1 + ctx->a * ctx->a * ctx->b * ctx->b);
 68:   VecRestoreArray(U, &u);
 69:   return 0;
 70: }

 72: int main(int argc, char **argv)
 73: {
 74:   TS           ts; /* ODE integrator */
 75:   Vec          U;  /* solution will be stored here */
 76:   Vec          Utrue;
 77:   PetscMPIInt  size;
 78:   AppCtx       ctx;
 79:   PetscScalar *u;
 80:   IS           iss;
 81:   IS           isf;
 82:   PetscInt    *indicess;
 83:   PetscInt    *indicesf;
 84:   PetscInt     n = 2;
 85:   PetscReal    error, tt;

 87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88:      Initialize program
 89:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 91:   PetscInitialize(&argc, &argv, (char *)0, help);
 92:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:     Create index for slow part and fast part
 97:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 98:   PetscMalloc1(1, &indicess);
 99:   indicess[0] = 0;
100:   PetscMalloc1(1, &indicesf);
101:   indicesf[0] = 1;
102:   ISCreateGeneral(PETSC_COMM_SELF, 1, indicess, PETSC_COPY_VALUES, &iss);
103:   ISCreateGeneral(PETSC_COMM_SELF, 1, indicesf, PETSC_COPY_VALUES, &isf);

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:     Create necessary vector
107:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108:   VecCreate(PETSC_COMM_WORLD, &U);
109:   VecSetSizes(U, n, PETSC_DETERMINE);
110:   VecSetFromOptions(U);
111:   VecDuplicate(U, &Utrue);
112:   VecCopy(U, Utrue);

114:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115:     Set runtime options
116:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ODE options", "");
118:   {
119:     ctx.a = 2.0;
120:     ctx.b = 25.0;
121:     PetscOptionsReal("-a", "", "", ctx.a, &ctx.a, NULL);
122:     PetscOptionsReal("-b", "", "", ctx.b, &ctx.b, NULL);
123:     ctx.Tf = 2;
124:     ctx.dt = 0.01;
125:     PetscOptionsReal("-Tf", "", "", ctx.Tf, &ctx.Tf, NULL);
126:     PetscOptionsReal("-dt", "", "", ctx.dt, &ctx.dt, NULL);
127:   }
128:   PetscOptionsEnd();

130:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131:      Initialize the solution
132:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133:   VecGetArray(U, &u);
134:   u[0] = 1.0;
135:   u[1] = ctx.a / (1 + ctx.a * ctx.a * ctx.b * ctx.b);
136:   VecRestoreArray(U, &u);

138:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139:      Create timestepping solver context
140:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141:   TSCreate(PETSC_COMM_WORLD, &ts);
142:   TSSetType(ts, TSMPRK);

144:   TSSetRHSFunction(ts, NULL, (TSRHSFunction)RHSFunction, &ctx);
145:   TSRHSSplitSetIS(ts, "slow", iss);
146:   TSRHSSplitSetIS(ts, "fast", isf);
147:   TSRHSSplitSetRHSFunction(ts, "slow", NULL, (TSRHSFunction)RHSFunctionslow, &ctx);
148:   TSRHSSplitSetRHSFunction(ts, "fast", NULL, (TSRHSFunction)RHSFunctionfast, &ctx);

150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:      Set initial conditions
152:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153:   TSSetSolution(ts, U);

155:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156:      Set solver options
157:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158:   TSSetMaxTime(ts, ctx.Tf);
159:   TSSetTimeStep(ts, ctx.dt);
160:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
161:   TSSetFromOptions(ts);

163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164:      Solve linear system
165:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166:   TSSolve(ts, U);
167:   VecView(U, PETSC_VIEWER_STDOUT_WORLD);

169:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170:      Check the error of the Petsc solution
171:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
172:   TSGetTime(ts, &tt);
173:   sol_true(tt, Utrue, &ctx);
174:   VecAXPY(Utrue, -1.0, U);
175:   VecNorm(Utrue, NORM_2, &error);

177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178:      Print norm2 error
179:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180:   PetscPrintf(PETSC_COMM_WORLD, "l2 error norm = %g\n", (double)error);

182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
184:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185:   VecDestroy(&U);
186:   VecDestroy(&Utrue);
187:   TSDestroy(&ts);
188:   ISDestroy(&iss);
189:   ISDestroy(&isf);
190:   PetscFree(indicess);
191:   PetscFree(indicesf);
192:   PetscFinalize();
193:   return 0;
194: }

196: /*TEST
197:     build:
198:       requires: !complex

200:     test:

202: TEST*/