Actual source code: ex24.c

  1: static char help[] = "Pseudotransient continuation to solve a many-variable system that comes from the 2 variable Rosenbrock function + trivial.\n\n";

  3: #include <petscts.h>

  5: static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
  6: static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
  7: static PetscErrorCode MonitorObjective(TS, PetscInt, PetscReal, Vec, void *);

  9: typedef struct {
 10:   PetscInt  n;
 11:   PetscBool monitor_short;
 12: } Ctx;

 14: int main(int argc, char **argv)
 15: {
 16:   TS           ts; /* time integration context */
 17:   Vec          X;  /* solution, residual vectors */
 18:   Mat          J;  /* Jacobian matrix */
 19:   PetscScalar *x;
 20:   PetscReal    ftime;
 21:   PetscInt     i, steps, nits, lits;
 22:   PetscBool    view_final;
 23:   Ctx          ctx;

 26:   PetscInitialize(&argc, &argv, (char *)0, help);
 27:   ctx.n = 3;
 28:   PetscOptionsGetInt(NULL, NULL, "-n", &ctx.n, NULL);

 31:   view_final = PETSC_FALSE;
 32:   PetscOptionsGetBool(NULL, NULL, "-view_final", &view_final, NULL);

 34:   ctx.monitor_short = PETSC_FALSE;
 35:   PetscOptionsGetBool(NULL, NULL, "-monitor_short", &ctx.monitor_short, NULL);

 37:   /*
 38:      Create Jacobian matrix data structure and state vector
 39:   */
 40:   MatCreate(PETSC_COMM_WORLD, &J);
 41:   MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, ctx.n, ctx.n);
 42:   MatSetFromOptions(J);
 43:   MatSetUp(J);
 44:   MatCreateVecs(J, &X, NULL);

 46:   /* Create time integration context */
 47:   TSCreate(PETSC_COMM_WORLD, &ts);
 48:   TSSetType(ts, TSPSEUDO);
 49:   TSSetIFunction(ts, NULL, FormIFunction, &ctx);
 50:   TSSetIJacobian(ts, J, J, FormIJacobian, &ctx);
 51:   TSSetMaxSteps(ts, 1000);
 52:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
 53:   TSSetTimeStep(ts, 1e-3);
 54:   TSMonitorSet(ts, MonitorObjective, &ctx, NULL);

 56:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 57:      Customize time integrator; set runtime options
 58:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 59:   TSSetFromOptions(ts);

 61:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 62:      Evaluate initial guess; then solve nonlinear system
 63:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 64:   VecSet(X, 0.0);
 65:   VecGetArray(X, &x);
 66: #if 1
 67:   x[0] = 5.;
 68:   x[1] = -5.;
 69:   for (i = 2; i < ctx.n; i++) x[i] = 5.;
 70: #else
 71:   x[0] = 1.0;
 72:   x[1] = 15.0;
 73:   for (i = 2; i < ctx.n; i++) x[i] = 10.0;
 74: #endif
 75:   VecRestoreArray(X, &x);

 77:   TSSolve(ts, X);
 78:   TSGetSolveTime(ts, &ftime);
 79:   TSGetStepNumber(ts, &steps);
 80:   TSGetSNESIterations(ts, &nits);
 81:   TSGetKSPIterations(ts, &lits);
 82:   PetscPrintf(PETSC_COMM_WORLD, "Time integrator took (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ") iterations to reach final time %g\n", steps, nits, lits, (double)ftime);
 83:   if (view_final) VecView(X, PETSC_VIEWER_STDOUT_WORLD);

 85:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 86:      Free work space.  All PETSc objects should be destroyed when they
 87:      are no longer needed.
 88:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 90:   VecDestroy(&X);
 91:   MatDestroy(&J);
 92:   TSDestroy(&ts);
 93:   PetscFinalize();
 94:   return 0;
 95: }

 97: static PetscErrorCode MonitorObjective(TS ts, PetscInt step, PetscReal t, Vec X, void *ictx)
 98: {
 99:   Ctx               *ctx = (Ctx *)ictx;
100:   const PetscScalar *x;
101:   PetscScalar        f;
102:   PetscReal          dt, gnorm;
103:   PetscInt           i, snesit, linit;
104:   SNES               snes;
105:   Vec                Xdot, F;

108:   /* Compute objective functional */
109:   VecGetArrayRead(X, &x);
110:   f = 0;
111:   for (i = 0; i < ctx->n - 1; i++) f += PetscSqr(1. - x[i]) + 100. * PetscSqr(x[i + 1] - PetscSqr(x[i]));
112:   VecRestoreArrayRead(X, &x);

114:   /* Compute norm of gradient */
115:   VecDuplicate(X, &Xdot);
116:   VecDuplicate(X, &F);
117:   VecZeroEntries(Xdot);
118:   FormIFunction(ts, t, X, Xdot, F, ictx);
119:   VecNorm(F, NORM_2, &gnorm);
120:   VecDestroy(&Xdot);
121:   VecDestroy(&F);

123:   TSGetTimeStep(ts, &dt);
124:   TSGetSNES(ts, &snes);
125:   SNESGetIterationNumber(snes, &snesit);
126:   SNESGetLinearSolveIterations(snes, &linit);
127:   PetscPrintf(PETSC_COMM_WORLD, (ctx->monitor_short ? "%3" PetscInt_FMT " t=%10.1e  dt=%10.1e  f=%10.1e  df=%10.1e  it=(%2" PetscInt_FMT ",%3" PetscInt_FMT ")\n" : "%3" PetscInt_FMT " t=%10.4e  dt=%10.4e  f=%10.4e  df=%10.4e  it=(%2" PetscInt_FMT ",%3" PetscInt_FMT ")\n"), step, (double)t, (double)dt, (double)PetscRealPart(f), (double)gnorm, snesit, linit);
128:   return 0;
129: }

131: /* ------------------------------------------------------------------- */
132: /*
133:    FormIFunction - Evaluates nonlinear function, F(X,Xdot) = Xdot + grad(objective(X))

135:    Input Parameters:
136: +  ts   - the TS context
137: .  t - time
138: .  X    - input vector
139: .  Xdot - time derivative
140: -  ctx  - optional user-defined context

142:    Output Parameter:
143: .  F - function vector
144:  */
145: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ictx)
146: {
147:   const PetscScalar *x;
148:   PetscScalar       *f;
149:   PetscInt           i;
150:   Ctx               *ctx = (Ctx *)ictx;

153:   /*
154:     Get pointers to vector data.
155:     - For default PETSc vectors, VecGetArray() returns a pointer to
156:     the data array.  Otherwise, the routine is implementation dependent.
157:     - You MUST call VecRestoreArray() when you no longer need access to
158:     the array.
159:   */
160:   VecGetArrayRead(X, &x);
161:   VecZeroEntries(F);
162:   VecGetArray(F, &f);

164:   /* Compute gradient of objective */
165:   for (i = 0; i < ctx->n - 1; i++) {
166:     PetscScalar a, a0, a1;
167:     a  = x[i + 1] - PetscSqr(x[i]);
168:     a0 = -2. * x[i];
169:     a1 = 1.;
170:     f[i] += -2. * (1. - x[i]) + 200. * a * a0;
171:     f[i + 1] += 200. * a * a1;
172:   }
173:   /* Restore vectors */
174:   VecRestoreArrayRead(X, &x);
175:   VecRestoreArray(F, &f);
176:   VecAXPY(F, 1.0, Xdot);
177:   return 0;
178: }
179: /* ------------------------------------------------------------------- */
180: /*
181:    FormIJacobian - Evaluates Jacobian matrix.

183:    Input Parameters:
184: +  ts - the TS context
185: .  t - pseudo-time
186: .  X - input vector
187: .  Xdot - time derivative
188: .  shift - multiplier for mass matrix
189: .  dummy - user-defined context

191:    Output Parameters:
192: .  J - Jacobian matrix
193: .  B - optionally different preconditioning matrix
194: .  flag - flag indicating matrix structure
195: */
196: static PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal shift, Mat J, Mat B, void *ictx)
197: {
198:   const PetscScalar *x;
199:   PetscInt           i;
200:   Ctx               *ctx = (Ctx *)ictx;

203:   MatZeroEntries(B);
204:   /*
205:      Get pointer to vector data
206:   */
207:   VecGetArrayRead(X, &x);

209:   /*
210:      Compute Jacobian entries and insert into matrix.
211:   */
212:   for (i = 0; i < ctx->n - 1; i++) {
213:     PetscInt    rowcol[2];
214:     PetscScalar v[2][2], a, a0, a1, a00, a01, a10, a11;
215:     rowcol[0] = i;
216:     rowcol[1] = i + 1;
217:     a         = x[i + 1] - PetscSqr(x[i]);
218:     a0        = -2. * x[i];
219:     a00       = -2.;
220:     a01       = 0.;
221:     a1        = 1.;
222:     a10       = 0.;
223:     a11       = 0.;
224:     v[0][0]   = 2. + 200. * (a * a00 + a0 * a0);
225:     v[0][1]   = 200. * (a * a01 + a1 * a0);
226:     v[1][0]   = 200. * (a * a10 + a0 * a1);
227:     v[1][1]   = 200. * (a * a11 + a1 * a1);
228:     MatSetValues(B, 2, rowcol, 2, rowcol, &v[0][0], ADD_VALUES);
229:   }
230:   for (i = 0; i < ctx->n; i++) MatSetValue(B, i, i, (PetscScalar)shift, ADD_VALUES);

232:   VecRestoreArrayRead(X, &x);

234:   /*
235:      Assemble matrix
236:   */
237:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
238:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
239:   if (J != B) {
240:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
241:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
242:   }
243:   return 0;
244: }

246: /*TEST

248:     test:
249:       requires: !single

251:     test:
252:       args: -pc_type lu -ts_dt 1e-5 -ts_max_time 1e5 -n 50 -monitor_short -snes_max_it 5 -snes_type newtonls -ts_max_snes_failures -1
253:       requires: !single
254:       suffix: 2

256: TEST*/