Actual source code: ex24.c

petsc-3.4.4 2014-03-13
  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*,MatStructure*,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;

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

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

 33:   view_final = PETSC_FALSE;
 34:   PetscOptionsGetBool(NULL,"-view_final",&view_final,NULL);

 36:   ctx.monitor_short = PETSC_FALSE;
 37:   PetscOptionsGetBool(NULL,"-monitor_short",&ctx.monitor_short,NULL);

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

 48:   /* Create time integration context */
 49:   TSCreate(PETSC_COMM_WORLD,&ts);
 50:   TSSetType(ts,TSPSEUDO);
 51:   TSSetIFunction(ts,NULL,FormIFunction,&ctx);
 52:   TSSetIJacobian(ts,J,J,FormIJacobian,&ctx);
 53:   TSSetDuration(ts,1000,1e14);
 54:   TSSetInitialTimeStep(ts,0.0,1e-3);
 55:   TSMonitorSet(ts,MonitorObjective,&ctx,NULL);

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

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

 78:   TSSolve(ts,X);
 79:   TSGetSolveTime(ts,&ftime);
 80:   TSGetTimeStepNumber(ts,&steps);
 81:   TSGetSNESIterations(ts,&nits);
 82:   TSGetKSPIterations(ts,&lits);
 83:   PetscPrintf(PETSC_COMM_WORLD,"Time integrator took (%D,%D,%D) iterations to reach final time %G\n",steps,nits,lits,ftime);
 84:   if (view_final) {
 85:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
 86:   }

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

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

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

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

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

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


142: /* ------------------------------------------------------------------- */
145: /*
146:    FormIFunction - Evaluates nonlinear function, F(X,Xdot) = Xdot + grad(objective(X))

148:    Input Parameters:
149: +  ts   - the TS context
150: .  t - time
151: .  X    - input vector
152: .  Xdot - time derivative
153: -  ctx  - optional user-defined context

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

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

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

199:    Input Parameters:
200: +  ts - the TS context
201: .  t - pseudo-time
202: .  X - input vector
203: .  Xdot - time derivative
204: .  shift - multiplier for mass matrix
205: .  dummy - user-defined context

207:    Output Parameters:
208: .  J - Jacobian matrix
209: .  B - optionally different preconditioning matrix
210: .  flag - flag indicating matrix structure
211: */
212: static PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat *J,Mat *B,MatStructure *flag,void *ictx)
213: {
214:   const PetscScalar *x;
215:   PetscErrorCode    ierr;
216:   PetscInt          i;
217:   Ctx               *ctx = (Ctx*)ictx;

220:   MatZeroEntries(*B);
221:   /*
222:      Get pointer to vector data
223:   */
224:   VecGetArrayRead(X,&x);

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

251:   VecRestoreArrayRead(X,&x);

253:   /*
254:      Assemble matrix
255:   */
256:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
257:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
258:   if (*J != *B) {
259:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
260:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
261:   }
262:   return(0);
263: }