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 */
 20:   PetscScalar    *x;
 21:   PetscReal      ftime;
 22:   PetscInt       i,steps,nits,lits;
 23:   PetscBool      view_final;
 24:   Ctx            ctx;

 26:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 27:   ctx.n = 3;
 28:   PetscOptionsGetInt(NULL,NULL,"-n",&ctx.n,NULL);
 29:   if (ctx.n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"The dimension specified with -n must be at least 2");

 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 (%D,%D,%D) iterations to reach final time %g\n",steps,nits,lits,(double)ftime);
 83:   if (view_final) {
 84:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
 85:   }

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

 92:   VecDestroy(&X);
 93:   MatDestroy(&J);
 94:   TSDestroy(&ts);
 95:   PetscFinalize();
 96:   return ierr;
 97: }

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

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

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

126:   TSGetTimeStep(ts,&dt);
127:   TSGetSNES(ts,&snes);
128:   SNESGetIterationNumber(snes,&snesit);
129:   SNESGetLinearSolveIterations(snes,&linit);
130:   PetscPrintf(PETSC_COMM_WORLD,
131:                      (ctx->monitor_short
132:                       ? "%3D t=%10.1e  dt=%10.1e  f=%10.1e  df=%10.1e  it=(%2D,%3D)\n"
133:                       : "%3D t=%10.4e  dt=%10.4e  f=%10.4e  df=%10.4e  it=(%2D,%3D)\n"),
134:                      step,(double)t,(double)dt,(double)PetscRealPart(f),(double)gnorm,snesit,linit);
135:   return(0);
136: }


139: /* ------------------------------------------------------------------- */
140: /*
141:    FormIFunction - Evaluates nonlinear function, F(X,Xdot) = Xdot + grad(objective(X))

143:    Input Parameters:
144: +  ts   - the TS context
145: .  t - time
146: .  X    - input vector
147: .  Xdot - time derivative
148: -  ctx  - optional user-defined context

150:    Output Parameter:
151: .  F - function vector
152:  */
153: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ictx)
154: {
155:   PetscErrorCode    ierr;
156:   const PetscScalar *x;
157:   PetscScalar       *f;
158:   PetscInt          i;
159:   Ctx               *ctx = (Ctx*)ictx;

162:   /*
163:     Get pointers to vector data.
164:     - For default PETSc vectors, VecGetArray() returns a pointer to
165:     the data array.  Otherwise, the routine is implementation dependent.
166:     - You MUST call VecRestoreArray() when you no longer need access to
167:     the array.
168:   */
169:   VecGetArrayRead(X,&x);
170:   VecZeroEntries(F);
171:   VecGetArray(F,&f);

173:   /* Compute gradient of objective */
174:   for (i=0; i<ctx->n-1; i++) {
175:     PetscScalar a,a0,a1;
176:     a       = x[i+1] - PetscSqr(x[i]);
177:     a0      = -2.*x[i];
178:     a1      = 1.;
179:     f[i]   += -2.*(1. - x[i]) + 200.*a*a0;
180:     f[i+1] += 200.*a*a1;
181:   }
182:   /* Restore vectors */
183:   VecRestoreArrayRead(X,&x);
184:   VecRestoreArray(F,&f);
185:   VecAXPY(F,1.0,Xdot);
186:   return(0);
187: }
188: /* ------------------------------------------------------------------- */
189: /*
190:    FormIJacobian - Evaluates Jacobian matrix.

192:    Input Parameters:
193: +  ts - the TS context
194: .  t - pseudo-time
195: .  X - input vector
196: .  Xdot - time derivative
197: .  shift - multiplier for mass matrix
198: .  dummy - user-defined context

200:    Output Parameters:
201: .  J - Jacobian matrix
202: .  B - optionally different preconditioning matrix
203: .  flag - flag indicating matrix structure
204: */
205: static PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat J,Mat B,void *ictx)
206: {
207:   const PetscScalar *x;
208:   PetscErrorCode    ierr;
209:   PetscInt          i;
210:   Ctx               *ctx = (Ctx*)ictx;

213:   MatZeroEntries(B);
214:   /*
215:      Get pointer to vector data
216:   */
217:   VecGetArrayRead(X,&x);

219:   /*
220:      Compute Jacobian entries and insert into matrix.
221:   */
222:   for (i=0; i<ctx->n-1; i++) {
223:     PetscInt    rowcol[2];
224:     PetscScalar v[2][2],a,a0,a1,a00,a01,a10,a11;
225:     rowcol[0] = i;
226:     rowcol[1] = i+1;
227:     a         = x[i+1] - PetscSqr(x[i]);
228:     a0        = -2.*x[i];
229:     a00       = -2.;
230:     a01       = 0.;
231:     a1        = 1.;
232:     a10       = 0.;
233:     a11       = 0.;
234:     v[0][0]   = 2. + 200.*(a*a00 + a0*a0);
235:     v[0][1]   = 200.*(a*a01 + a1*a0);
236:     v[1][0]   = 200.*(a*a10 + a0*a1);
237:     v[1][1]   = 200.*(a*a11 + a1*a1);
238:     MatSetValues(B,2,rowcol,2,rowcol,&v[0][0],ADD_VALUES);
239:   }
240:   for (i=0; i<ctx->n; i++) {
241:     MatSetValue(B,i,i,(PetscScalar)shift,ADD_VALUES);
242:   }

244:   VecRestoreArrayRead(X,&x);

246:   /*
247:      Assemble matrix
248:   */
249:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
250:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
251:   if (J != B) {
252:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
253:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
254:   }
255:   return(0);
256: }

258: /*TEST

260:     test:
261:       requires: !single

263:     test:
264:       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
265:       requires: !single
266:       suffix: 2

268: TEST*/