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;

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

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

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

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

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

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

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

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

 81:   /* free the memories */

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

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

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

 97:   PetscFunctionBeginUser;
 98:   PetscCall(VecGetArrayRead(x, &inptr));
 99:   PetscCall(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:   PetscCall(VecRestoreArrayRead(x, &inptr));
106:   PetscCall(VecRestoreArrayWrite(y, &outptr));
107:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscFunctionBeginUser;
117:   /* determine starting point of each processor */
118:   PetscCall(VecGetOwnershipRange(global, &mybase, &myend));
119:   PetscCall(VecGetLocalSize(global, &locsize));

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

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

127:   PetscCall(VecRestoreArrayWrite(global, &localptr));
128:   PetscFunctionReturn(PETSC_SUCCESS);
129: }

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

139:   PetscFunctionBeginUser;
140:   /* Get the size of the vector */
141:   PetscCall(VecGetSize(global, &n));

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

147:   /* Create local sequential vectors */
148:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_vec));

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

157:   PetscCall(VecGetArrayRead(tmp_vec, &tmp));
158:   PetscCall(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])));
159:   PetscCall(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))));
160:   PetscCall(VecRestoreArrayRead(tmp_vec, &tmp));
161:   PetscCall(VecScatterDestroy(&scatter));
162:   PetscCall(ISDestroy(&from));
163:   PetscCall(ISDestroy(&to));
164:   PetscCall(PetscFree(idx));
165:   PetscCall(VecDestroy(&tmp_vec));
166:   PetscFunctionReturn(PETSC_SUCCESS);
167: }

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

178:   PetscFunctionBeginUser;
179:   /* Get the length of parallel vector */
180:   PetscCall(VecGetSize(globalin, &n));

182:   /* Set the index sets */
183:   PetscCall(PetscMalloc1(n, &idx));
184:   for (i = 0; i < n; i++) idx[i] = i;

186:   /* Create local sequential vectors */
187:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_in));
188:   PetscCall(VecDuplicate(tmp_in, &tmp_out));

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

198:   /*Extract income array */
199:   PetscCall(VecGetArrayRead(tmp_in, &inptr));

201:   /* Extract outcome array*/
202:   PetscCall(VecGetArrayWrite(tmp_out, &outptr));

204:   outptr[0] = 2.0 * inptr[0] + inptr[1];
205:   outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
206:   outptr[2] = inptr[1] + 2.0 * inptr[2];

208:   PetscCall(VecRestoreArrayRead(tmp_in, &inptr));
209:   PetscCall(VecRestoreArrayWrite(tmp_out, &outptr));

211:   PetscCall(VecScatterCreate(tmp_out, from, globalout, to, &scatter));
212:   PetscCall(VecScatterBegin(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD));
213:   PetscCall(VecScatterEnd(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD));

215:   /* Destroy idx aand scatter */
216:   PetscCall(ISDestroy(&from));
217:   PetscCall(ISDestroy(&to));
218:   PetscCall(VecScatterDestroy(&scatter));
219:   PetscCall(VecDestroy(&tmp_in));
220:   PetscCall(VecDestroy(&tmp_out));
221:   PetscCall(PetscFree(idx));
222:   PetscFunctionReturn(PETSC_SUCCESS);
223: }

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

231:   PetscFunctionBeginUser;
232:   idx[0] = 0;
233:   idx[1] = 1;
234:   idx[2] = 2;
235:   PetscCall(VecGetArrayRead(x, &tmp));

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

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

249:   i    = 2;
250:   v[0] = 0.0;
251:   v[1] = 1.0;
252:   v[2] = 2.0;
253:   PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES));

255:   PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY));
256:   PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY));

258:   if (A != BB) {
259:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
260:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
261:   }
262:   PetscCall(VecRestoreArrayRead(x, &tmp));

264:   PetscFunctionReturn(PETSC_SUCCESS);
265: }

267: /*
268:       The exact solutions
269: */
270: PetscReal solx(PetscReal t)
271: {
272:   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));
273: }

275: PetscReal soly(PetscReal t)
276: {
277:   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);
278: }

280: PetscReal solz(PetscReal t)
281: {
282:   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));
283: }

285: /*TEST

287:     test:
288:       suffix: euler
289:       args: -ts_type euler
290:       requires: !single

292:     test:
293:       suffix: beuler
294:       args:   -ts_type beuler
295:       requires: !single

297: TEST*/