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*/