Actual source code: ex43.c

petsc-3.4.4 2014-03-13
  2: static char help[] = "Newton's method to solve a many-variable system that comes from the 2 variable Rosenbrock function + trivial.\n\n";

  4: /*

  6: ./ex43 -snes_monitor_range -snes_max_it 1000 -snes_rtol 1.e-14 -n 10 -snes_converged_reason -sub_snes_monito -sub_snes_mf -sub_snes_converged_reason -sub_snes_rtol 1.e-10 -sub_snes_max_it 1000 -sub_snes_monitor

  8:   Accelerates Newton's method by solving a small problem defined by those elements with large residual plus one level of overlap

 10:   This is a toy code for playing around

 12:   Counts residual entries as small if they are less then .2 times the maximum
 13:   Decides to solve a reduced problem if the number of large entries is less than 20 percent of all entries (and this has been true for criteria_reduce iterations)
 14: */
 15: #include "ex43-44.h"


 18: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 19: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);

 21: typedef struct {
 22:   PetscInt n,p;
 23: } Ctx;

 27: int main(int argc,char **argv)
 28: {
 29:   SNES                snes;         /* nonlinear solver context */
 30:   Vec                 x,r;          /* solution, residual vectors */
 31:   Mat                 J;            /* Jacobian matrix */
 32:   PetscErrorCode      ierr;
 33:   PetscScalar         *xx;
 34:   PetscInt            i,max_snes_solves = 20,snes_steps_per_solve = 2,criteria_reduce = 1;
 35:   Ctx                 ctx;
 36:   SNESConvergedReason reason;

 38:   PetscInitialize(&argc,&argv,(char*)0,help);
 39:   ctx.n = 0;
 40:   PetscOptionsGetInt(NULL,"-n",&ctx.n,NULL);
 41:   ctx.p = 0;
 42:   PetscOptionsGetInt(NULL,"-p",&ctx.p,NULL);
 43:   PetscOptionsGetInt(NULL,"-max_snes_solves",&max_snes_solves,NULL);
 44:   PetscOptionsGetInt(NULL,"-snes_steps_per_solve",&snes_steps_per_solve,NULL);
 45:   PetscOptionsGetInt(NULL,"-criteria_reduce",&criteria_reduce,NULL);

 47:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 48:      Create nonlinear solver context
 49:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 50:   SNESCreate(PETSC_COMM_WORLD,&snes);

 52:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53:      Create matrix and vector data structures; set corresponding routines
 54:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 55:   /*
 56:      Create vectors for solution and nonlinear function
 57:   */
 58:   VecCreate(PETSC_COMM_WORLD,&x);
 59:   VecSetSizes(x,PETSC_DECIDE,2+ctx.n+ctx.p);
 60:   VecSetFromOptions(x);
 61:   VecDuplicate(x,&r);

 63:   /*
 64:      Create Jacobian matrix data structure
 65:   */
 66:   MatCreate(PETSC_COMM_WORLD,&J);
 67:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2+ctx.p+ctx.n,2+ctx.p+ctx.n);
 68:   MatSetFromOptions(J);
 69:   MatSetUp(J);

 71:   /*
 72:      Set function evaluation routine and vector.
 73:   */
 74:   SNESSetFunction(snes,r,FormFunction1,(void*)&ctx);

 76:   /*
 77:      Set Jacobian matrix data structure and Jacobian evaluation routine
 78:   */
 79:   SNESSetJacobian(snes,J,J,FormJacobian1,(void*)&ctx);

 81:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 82:      Customize nonlinear solver; set runtime options
 83:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 84:   SNESSetFromOptions(snes);

 86:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 87:      Evaluate initial guess; then solve nonlinear system
 88:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 89:   VecSet(x,0.0);
 90:   VecGetArray(x,&xx);
 91:   xx[0] = -1.2;
 92:   for (i=1; i<ctx.p+2; i++) xx[i] = 1.0;
 93:   VecRestoreArray(x,&xx);

 95:   /*
 96:      Note: The user should initialize the vector, x, with the initial guess
 97:      for the nonlinear solver prior to calling SNESSolve().  In particular,
 98:      to employ an initial guess of zero, the user should explicitly set
 99:      this vector to zero by calling VecSet().
100:   */

102:   SNESMonitorSet(snes,MonitorRange,0,0);
103:   SNESSetTolerances(snes,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,snes_steps_per_solve,PETSC_DEFAULT);
104:   for (i=0; i<max_snes_solves; i++) {
105:     SNESSolve(snes,NULL,x);
106:     SNESGetConvergedReason(snes,&reason);
107:     if (reason && reason != SNES_DIVERGED_MAX_IT) break;
108:     if (CountGood > criteria_reduce) {
109:       SolveSubproblem(snes);
110:       CountGood = 0;
111:     }
112:   }

114:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115:      Free work space.  All PETSc objects should be destroyed when they
116:      are no longer needed.
117:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

119:   VecDestroy(&x); VecDestroy(&r);
120:   MatDestroy(&J); SNESDestroy(&snes);
121:   PetscFinalize();
122:   return 0;
123: }
124: /* ------------------------------------------------------------------- */
127: /*
128:    FormFunction1 - Evaluates nonlinear function, F(x).

130:    Input Parameters:
131: .  snes - the SNES context
132: .  x    - input vector
133: .  ctx  - optional user-defined context

135:    Output Parameter:
136: .  f - function vector
137:  */
138: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ictx)
139: {
141:   PetscScalar    *xx,*ff;
142:   PetscInt       i;
143:   Ctx            *ctx = (Ctx*)ictx;

145:   /*
146:     Get pointers to vector data.
147:     - For default PETSc vectors, VecGetArray() returns a pointer to
148:     the data array.  Otherwise, the routine is implementation dependent.
149:     - You MUST call VecRestoreArray() when you no longer need access to
150:     the array.
151:   */
152:   VecGetArray(x,&xx);
153:   VecGetArray(f,&ff);

155:   /* Compute function */
156:   ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1];
157:   for (i=1; i<1+ctx->p; i++) {
158:     ff[i] = -2.0 + 2.0*xx[i] + 400.0*xx[i]*xx[i]*xx[i] - 400.0*xx[i]*xx[i+1] + 200.0*(xx[i] - xx[i-1]*xx[i-1]);
159:   }
160:   ff[ctx->p+1] = -200.0*xx[ctx->p]*xx[ctx->p] + 200.0*xx[ctx->p+1];
161:   for (i=ctx->p+2; i<2+ctx->p+ctx->n; i++) {
162:     ff[i] = xx[i] - xx[0] + .7*xx[1] - .2*xx[i-1]*xx[i-1];
163:   }

165:   /* Restore vectors */
166:   VecRestoreArray(x,&xx);
167:   VecRestoreArray(f,&ff);
168:   return 0;
169: }
170: /* ------------------------------------------------------------------- */
173: /*
174:    FormJacobian1 - Evaluates Jacobian matrix.

176:    Input Parameters:
177: .  snes - the SNES context
178: .  x - input vector
179: .  dummy - optional user-defined context (not used here)

181:    Output Parameters:
182: .  jac - Jacobian matrix
183: .  B - optionally different preconditioning matrix
184: .  flag - flag indicating matrix structure
185: */
186: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *ictx)
187: {
188:   PetscScalar    *xx;
190:   PetscInt       i;
191:   Ctx            *ctx = (Ctx*)ictx;

193:   MatZeroEntries(*B);
194:   /*
195:      Get pointer to vector data
196:   */
197:   VecGetArray(x,&xx);

199:   /*
200:      Compute Jacobian entries and insert into matrix.
201:       - Since this is such a small problem, we set all entries for
202:         the matrix at once.
203:   */
204:   MatSetValue(*B,0,0, 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1],ADD_VALUES);
205:   MatSetValue(*B,0,1,-400.0*xx[0],ADD_VALUES);

207:   for (i=1; i<ctx->p+1; i++) {
208:     MatSetValue(*B,i,i-1, -400.0*xx[i-1],ADD_VALUES);
209:     MatSetValue(*B,i,i, 2.0 + 1200.0*xx[i]*xx[i] - 400.0*xx[i+1] + 200.0,ADD_VALUES);
210:     MatSetValue(*B,i,i+1,-400.0*xx[i],ADD_VALUES);
211:   }

213:   MatSetValue(*B,ctx->p+1,ctx->p, -400.0*xx[ctx->p],ADD_VALUES);
214:   MatSetValue(*B,ctx->p+1,ctx->p+1,200,ADD_VALUES);

216:   *flag = SAME_NONZERO_PATTERN;

218:   for (i=ctx->p+2; i<2+ctx->p+ctx->n; i++) {
219:     MatSetValue(*B,i,i,1.0,ADD_VALUES);
220:     MatSetValue(*B,i,0,-1.0,ADD_VALUES);
221:     MatSetValue(*B,i,1,.7,ADD_VALUES);
222:     MatSetValue(*B,i,i-1,-.4*xx[i-1],ADD_VALUES);
223:   }
224:   /*
225:      Restore vector
226:   */
227:   VecRestoreArray(x,&xx);

229:   /*
230:      Assemble matrix
231:   */
232:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
233:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
234:   if (*jac != *B) {
235:     MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
236:     MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
237:   }
238:   return 0;
239: }