Actual source code: ex2.c

  1: /*
  2:        Formatted test for TS routines.

  4:           Solves U_t=F(t,u)
  5:           Where:

  7:                   [2*u1+u2
  8:           F(t,u)= [u1+2*u2+u3
  9:                   [   u2+2*u3
 10:        We can compare the solutions from euler, beuler and SUNDIALS to
 11:        see what is the difference.

 13: */

 15: static char help[] = "Solves a linear ODE. \n\n";

 17: #include <petscts.h>
 18: #include <petscpc.h>

 20: extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
 21: extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
 22: extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
 23: extern PetscErrorCode Initial(Vec, void *);
 24: extern PetscErrorCode MyMatMult(Mat, Vec, Vec);

 26: extern PetscReal solx(PetscReal);
 27: extern PetscReal soly(PetscReal);
 28: extern PetscReal solz(PetscReal);

 30: int main(int argc, char **argv)
 31: {
 32:   PetscInt  time_steps = 100, steps;
 33:   Vec       global;
 34:   PetscReal dt, ftime;
 35:   TS        ts;
 36:   Mat       A = 0, S;

 39:   PetscInitialize(&argc, &argv, (char *)0, help);
 40:   PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL);

 42:   /* set initial conditions */
 43:   VecCreate(PETSC_COMM_WORLD, &global);
 44:   VecSetSizes(global, PETSC_DECIDE, 3);
 45:   VecSetFromOptions(global);
 46:   Initial(global, NULL);

 48:   /* make timestep context */
 49:   TSCreate(PETSC_COMM_WORLD, &ts);
 50:   TSSetProblemType(ts, TS_NONLINEAR);
 51:   TSMonitorSet(ts, Monitor, NULL, NULL);
 52:   dt = 0.001;

 54:   /*
 55:     The user provides the RHS and Jacobian
 56:   */
 57:   TSSetRHSFunction(ts, NULL, RHSFunction, NULL);
 58:   MatCreate(PETSC_COMM_WORLD, &A);
 59:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 3, 3);
 60:   MatSetFromOptions(A);
 61:   MatSetUp(A);
 62:   RHSJacobian(ts, 0.0, global, A, A, NULL);
 63:   TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL);

 65:   MatCreateShell(PETSC_COMM_WORLD, 3, 3, 3, 3, NULL, &S);
 66:   MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MyMatMult);
 67:   TSSetRHSJacobian(ts, S, A, RHSJacobian, NULL);

 69:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
 70:   TSSetFromOptions(ts);

 72:   TSSetTimeStep(ts, dt);
 73:   TSSetMaxSteps(ts, time_steps);
 74:   TSSetMaxTime(ts, 1);
 75:   TSSetSolution(ts, global);

 77:   TSSolve(ts, global);
 78:   TSGetSolveTime(ts, &ftime);
 79:   TSGetStepNumber(ts, &steps);

 81:   /* free the memories */

 83:   TSDestroy(&ts);
 84:   VecDestroy(&global);
 85:   MatDestroy(&A);
 86:   MatDestroy(&S);

 88:   PetscFinalize();
 89:   return 0;
 90: }

 92: PetscErrorCode MyMatMult(Mat S, Vec x, Vec y)
 93: {
 94:   const PetscScalar *inptr;
 95:   PetscScalar       *outptr;

 98:   VecGetArrayRead(x, &inptr);
 99:   VecGetArrayWrite(y, &outptr);

101:   outptr[0] = 2.0 * inptr[0] + inptr[1];
102:   outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
103:   outptr[2] = inptr[1] + 2.0 * inptr[2];

105:   VecRestoreArrayRead(x, &inptr);
106:   VecRestoreArrayWrite(y, &outptr);
107:   return 0;
108: }

110: /* this test problem has initial values (1,1,1).                      */
111: PetscErrorCode Initial(Vec global, void *ctx)
112: {
113:   PetscScalar *localptr;
114:   PetscInt     i, mybase, myend, locsize;

116:   /* determine starting point of each processor */
117:   VecGetOwnershipRange(global, &mybase, &myend);
118:   VecGetLocalSize(global, &locsize);

120:   /* Initialize the array */
121:   VecGetArrayWrite(global, &localptr);
122:   for (i = 0; i < locsize; i++) localptr[i] = 1.0;

124:   if (mybase == 0) localptr[0] = 1.0;

126:   VecRestoreArrayWrite(global, &localptr);
127:   return 0;
128: }

130: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, void *ctx)
131: {
132:   VecScatter         scatter;
133:   IS                 from, to;
134:   PetscInt           i, n, *idx;
135:   Vec                tmp_vec;
136:   const PetscScalar *tmp;

138:   /* Get the size of the vector */
139:   VecGetSize(global, &n);

141:   /* Set the index sets */
142:   PetscMalloc1(n, &idx);
143:   for (i = 0; i < n; i++) idx[i] = i;

145:   /* Create local sequential vectors */
146:   VecCreateSeq(PETSC_COMM_SELF, n, &tmp_vec);

148:   /* Create scatter context */
149:   ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from);
150:   ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to);
151:   VecScatterCreate(global, from, tmp_vec, to, &scatter);
152:   VecScatterBegin(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD);
153:   VecScatterEnd(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD);

155:   VecGetArrayRead(tmp_vec, &tmp);
156:   PetscPrintf(PETSC_COMM_WORLD, "At t =%14.6e u = %14.6e  %14.6e  %14.6e \n", (double)time, (double)PetscRealPart(tmp[0]), (double)PetscRealPart(tmp[1]), (double)PetscRealPart(tmp[2]));
157:   PetscPrintf(PETSC_COMM_WORLD, "At t =%14.6e errors = %14.6e  %14.6e  %14.6e \n", (double)time, (double)PetscRealPart(tmp[0] - solx(time)), (double)PetscRealPart(tmp[1] - soly(time)), (double)PetscRealPart(tmp[2] - solz(time)));
158:   VecRestoreArrayRead(tmp_vec, &tmp);
159:   VecScatterDestroy(&scatter);
160:   ISDestroy(&from);
161:   ISDestroy(&to);
162:   PetscFree(idx);
163:   VecDestroy(&tmp_vec);
164:   return 0;
165: }

167: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)
168: {
169:   PetscScalar       *outptr;
170:   const PetscScalar *inptr;
171:   PetscInt           i, n, *idx;
172:   IS                 from, to;
173:   VecScatter         scatter;
174:   Vec                tmp_in, tmp_out;

176:   /* Get the length of parallel vector */
177:   VecGetSize(globalin, &n);

179:   /* Set the index sets */
180:   PetscMalloc1(n, &idx);
181:   for (i = 0; i < n; i++) idx[i] = i;

183:   /* Create local sequential vectors */
184:   VecCreateSeq(PETSC_COMM_SELF, n, &tmp_in);
185:   VecDuplicate(tmp_in, &tmp_out);

187:   /* Create scatter context */
188:   ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from);
189:   ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to);
190:   VecScatterCreate(globalin, from, tmp_in, to, &scatter);
191:   VecScatterBegin(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD);
192:   VecScatterEnd(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD);
193:   VecScatterDestroy(&scatter);

195:   /*Extract income array */
196:   VecGetArrayRead(tmp_in, &inptr);

198:   /* Extract outcome array*/
199:   VecGetArrayWrite(tmp_out, &outptr);

201:   outptr[0] = 2.0 * inptr[0] + inptr[1];
202:   outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
203:   outptr[2] = inptr[1] + 2.0 * inptr[2];

205:   VecRestoreArrayRead(tmp_in, &inptr);
206:   VecRestoreArrayWrite(tmp_out, &outptr);

208:   VecScatterCreate(tmp_out, from, globalout, to, &scatter);
209:   VecScatterBegin(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD);
210:   VecScatterEnd(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD);

212:   /* Destroy idx aand scatter */
213:   ISDestroy(&from);
214:   ISDestroy(&to);
215:   VecScatterDestroy(&scatter);
216:   VecDestroy(&tmp_in);
217:   VecDestroy(&tmp_out);
218:   PetscFree(idx);
219:   return 0;
220: }

222: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, void *ctx)
223: {
224:   PetscScalar        v[3];
225:   const PetscScalar *tmp;
226:   PetscInt           idx[3], i;

228:   idx[0] = 0;
229:   idx[1] = 1;
230:   idx[2] = 2;
231:   VecGetArrayRead(x, &tmp);

233:   i    = 0;
234:   v[0] = 2.0;
235:   v[1] = 1.0;
236:   v[2] = 0.0;
237:   MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES);

239:   i    = 1;
240:   v[0] = 1.0;
241:   v[1] = 2.0;
242:   v[2] = 1.0;
243:   MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES);

245:   i    = 2;
246:   v[0] = 0.0;
247:   v[1] = 1.0;
248:   v[2] = 2.0;
249:   MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES);

251:   MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY);
252:   MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY);

254:   if (A != BB) {
255:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
256:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
257:   }
258:   VecRestoreArrayRead(x, &tmp);

260:   return 0;
261: }

263: /*
264:       The exact solutions
265: */
266: PetscReal solx(PetscReal t)
267: {
268:   return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0)) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0));
269: }

271: PetscReal soly(PetscReal t)
272: {
273:   return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / PetscSqrtReal(2.0) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / PetscSqrtReal(2.0);
274: }

276: PetscReal solz(PetscReal t)
277: {
278:   return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0)) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0));
279: }

281: /*TEST

283:     test:
284:       suffix: euler
285:       args: -ts_type euler
286:       requires: !single

288:     test:
289:       suffix: beuler
290:       args:   -ts_type beuler
291:       requires: !single

293: TEST*/