Actual source code: ex3.c

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

  3: /*F

  5: \begin{eqnarray}
  6:                  ys' = -2.0*\frac{-1.0+ys^2.0-\cos(t)}{2.0*ys}+0.05*\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-\frac{\sin(t)}{2.0*ys}\\
  7:                  yf' = 0.05*\frac{-1.0+ys^2-\cos(t)}{2.0*ys}-\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-5.0*\frac{\sin(5.0*t)}{2.0*yf}\\
  8: \end{eqnarray}

 10: F*/

 12: #include <petscts.h>

 14: typedef struct {
 15:   PetscReal 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] = -2.0 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) + 0.05 * (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]) - PetscSinScalar(t) / (2.0 * u[0]);
 26:   f[1] = 0.05 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) - (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]) - 5.0 * PetscSinScalar(5.0 * t) / (2.0 * u[1]);
 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] = -2.0 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) + 0.05 * (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]) - PetscSinScalar(t) / (2.0 * u[0]);
 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] = 0.05 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) - (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]) - 5.0 * PetscSinScalar(5.0 * t) / (2.0 * u[1]);
 53:   VecRestoreArrayRead(U, &u);
 54:   VecRestoreArray(F, &f);
 55:   return 0;
 56: }

 58: static PetscErrorCode sol_true(PetscReal t, Vec U)
 59: {
 60:   PetscScalar *u;

 62:   VecGetArray(U, &u);
 63:   u[0] = PetscSqrtScalar(1.0 + PetscCosScalar(t));
 64:   u[1] = PetscSqrtScalar(2.0 + PetscCosScalar(5.0 * t));
 65:   VecRestoreArray(U, &u);
 66:   return 0;
 67: }

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

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85:      Initialize program
 86:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 88:   PetscInitialize(&argc, &argv, (char *)0, help);
 89:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

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

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

111:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112:     Set initial condition
113:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114:   VecGetArray(U, &u);
115:   u[0] = PetscSqrtScalar(2.0);
116:   u[1] = PetscSqrtScalar(3.0);
117:   VecRestoreArray(U, &u);

119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120:      Create timestepping solver context
121:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122:   TSCreate(PETSC_COMM_WORLD, &ts);
123:   TSSetType(ts, TSMPRK);

125:   TSSetRHSFunction(ts, NULL, (TSRHSFunction)RHSFunction, &ctx);
126:   TSRHSSplitSetIS(ts, "slow", iss);
127:   TSRHSSplitSetIS(ts, "fast", isf);
128:   TSRHSSplitSetRHSFunction(ts, "slow", NULL, (TSRHSFunction)RHSFunctionslow, &ctx);
129:   TSRHSSplitSetRHSFunction(ts, "fast", NULL, (TSRHSFunction)RHSFunctionfast, &ctx);

131:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132:      Set initial conditions
133:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134:   TSSetSolution(ts, U);

136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137:      Set solver options
138:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ODE options", "");
140:   {
141:     ctx.Tf = 0.3;
142:     ctx.dt = 0.01;
143:     PetscOptionsScalar("-Tf", "", "", ctx.Tf, &ctx.Tf, NULL);
144:     PetscOptionsScalar("-dt", "", "", ctx.dt, &ctx.dt, NULL);
145:   }
146:   PetscOptionsEnd();
147:   TSSetMaxTime(ts, ctx.Tf);
148:   TSSetTimeStep(ts, ctx.dt);
149:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
150:   TSSetFromOptions(ts);

152:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153:      Solve linear system
154:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155:   TSSolve(ts, U);
156:   VecView(U, PETSC_VIEWER_STDOUT_WORLD);

158:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159:      Check the error of the Petsc solution
160:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161:   TSGetTime(ts, &tt);
162:   sol_true(tt, Utrue);
163:   VecAXPY(Utrue, -1.0, U);
164:   VecNorm(Utrue, NORM_2, &error);

166:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167:      Print norm2 error
168:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169:   PetscPrintf(PETSC_COMM_WORLD, "l2 error norm: %g\n", (double)error);

171:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
173:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174:   VecDestroy(&U);
175:   TSDestroy(&ts);
176:   VecDestroy(&Utrue);
177:   ISDestroy(&iss);
178:   ISDestroy(&isf);
179:   PetscFree(indicess);
180:   PetscFree(indicesf);
181:   PetscFinalize();
182:   return 0;
183: }

185: /*TEST
186:     build:
187:       requires: !complex

189:     test:

191: TEST*/