Actual source code: ex43.c

petsc-3.3-p7 2013-05-11
```  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(PETSC_NULL,"-n",&ctx.n,PETSC_NULL);
41:   ctx.p = 0;
42:   PetscOptionsGetInt(PETSC_NULL,"-p",&ctx.p,PETSC_NULL);
43:   PetscOptionsGetInt(PETSC_NULL,"-max_snes_solves",&max_snes_solves,PETSC_NULL);
44:   PetscOptionsGetInt(PETSC_NULL,"-snes_steps_per_solve",&snes_steps_per_solve,PETSC_NULL);
45:   PetscOptionsGetInt(PETSC_NULL,"-criteria_reduce",&criteria_reduce,PETSC_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);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

217:   *flag = SAME_NONZERO_PATTERN;

219:   for (i=ctx->p+2; i<2+ctx->p+ctx->n; i++) {
224:   }
225:   /*
226:      Restore vector
227:   */
228:   VecRestoreArray(x,&xx);

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

```